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Numerical results for ground state and excited state properties (energies, double occupancies, 
and Matsubara-axis self energies) of the single-orbital Hubbard model on a two-dimensional square 
lattice are presented, in order to provide an assessment of our ability to compute accurate results in 
the thermodynamic limit. Many methods are employed, including auxiliary field quantum Monte 
Carlo, bare and bold-line diagrammatic Monte Carlo, method of dual fermions, density matrix 
embedding theory, density matrix renormalization group, dynamical cluster approximation, diffusion 
Monte Carlo within a fixed node approximation, unrestricted coupled cluster theory, and multi¬ 
reference projected Hartree-Fock. Comparison of results obtained by different methods allows for the 
identification of uncertainties and systematic errors. The importance of extrapolation to converged 
thermodynamic limit values is emphasized. Cases where agreement between different methods is 
obtained establish benchmark results that may be useful in the validation of new approaches and 
the improvement of existing methods. 


I. INTRODUCTION 

The “many-electron problem” of providing a useful and 
sufficiently accurate calculation of the properties of sys¬ 
tems of large numbers of interacting electrons is one of 
the grand scientific challenges of the present day. Im¬ 
proved solutions are needed both for the practical prob¬ 
lems of materials science and chemistry and for the basic 
science questions of determining the qualitative behav¬ 
iors of interacting quantum systems. 

While many problems of implementation arise, includ¬ 
ing calculation of the multiplicity of orbitals and inter¬ 
action matrix elements needed to characterize real mate¬ 
rials, the fundamental difficulties are that the dimension 
of the Hilbert space needed to describe an interacting 
electron system grows exponentially in the system size, 
so that direct diagonalization is not practical except for 
small systems, and that the minus sign associated with 
the Fermi statistics of electrons leads to exponentially 
slow convergence of straightforward Monte Carlo calcu¬ 


lations. It is generally accepted that a complete solu¬ 
tion to the many-electron problem cannot be obtained in 
polynomial time. 

The difficulties associated with obtaining a complete 
solution have motivated the development over the years 
of approximate methods, and comparison of the dif¬ 
ferent approximations remains a crucial open question. 
In this paper we address this issne in the context of 
the repulsive-interaction Hubbard modeP^^ defined on 
a two-dimensional square lattice. The Hubbard model 
is one of the simplest models of interacting fermions, 
but despite its simplicity exhibits a wide range of cor¬ 
related electron behavior including interaction-driven 
metal-insulator transitions, superconductivity, and mag¬ 
netism. The precise behavior depends delicately on pa¬ 
rameters, creating an interesting challenge for theory and 
computation. 

Exact solutions are ava ilab le for one-dimensionaP^ and 
infinite-dimensional cases.^^High temperature series ex¬ 
pansions provide numerically exact results, but only for 
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temperatures too h^h to be relevant for physically in¬ 
teresting situations^ In general dimensions at relevant 
temperatures, only approximate solutions are available. 
In some cases, these provide rigorous bounds bounds 
to the gr oun d state energies or other thermodynamic 
propertiesAnalytical perturbative methods can pro¬ 
vide information about the behavior at very small in¬ 
teraction strengtlPSHUl and at very large interaction 
stre ngth ( for the special case of nearly one electron per 
site) but outside of these limits obtaining results 
requires numerical methods.l^ Other techniques such as 
diagrammatic resummation are expected to work well in 
the weak coupling regime.^ The known numerical meth¬ 
ods are based on approximation schemes. Among the 
approximations employed are study of finite systems (ei¬ 
ther directly or via embedding constructions), use of vari¬ 
ational wave functions and evaluation of subsets of all 
possible Feynman diagrams. Controlling these approx¬ 
imations and assessing the remaining uncertainties is a 
challenging but essential task, requiring analysis of re¬ 
sults obtained from different methods. The past decade 
has seen the development of interesting new methods and 
substantial improvements in capabilities of previously de¬ 
veloped approaches, suggesting that the time is ripe for 
a careful comparison. 

In this paper we undertake this needed task. Our aim 
is to assess the state of our knowledge of the Hubbard 
model, identifying parameter regimes where reliable re¬ 
sults have been established and regimes where further 
work is needed. In regimes where reliable results have 
been established our results will serve as benchmarks to 
aid in the development and validation of new methods. 
This improved understanding of the Hubbard model will 
serve as a tool to analyze methods for solving the general 
many electron problem. 

We take the view that the only meaningful points of 
comparison between methods are results for the actual 
thermodynamic limit of the Hubbard model, with the 
uncertainties arising from the needed extrapolations (to 
infinite system size, to all diagrams, to infinite statis¬ 
tical precision in Monte Carlo calculations etc) quanti¬ 
fied. However, examination of results obtained at dif¬ 
ferent stages of the extrapolation sequence for a given 
method provide considerable information. Therefore, we 
present, when needed, both converged values and the in¬ 
termediate results from which these were obtained. 

The methods considered are auxiliary field quan¬ 
tum Monte Carlo (AFQMC),^2lHl3l bold-line 

diagrammatic Monte Carlo (DiagMC),l23H2ll the dual 
fermion met hod (D F),^^ density matrix embedding the¬ 
ory (DMET)) ^^ I ^ ^ I den sity matrix renormalization group 
theory (DMRG) cluster dynamical mean field the¬ 
ory in the dynamical cluster approximation (DCA),I^ 
diffusion Monte Carlo based on a fixed node approx¬ 
imation (FN),l22H3ni unrestricted coupled cluster theory 
including singles and doubles (UCCSD), and in certain 
cases, higher excitations, guj multi-reference pro¬ 
jected Hartree-Fock (MRPHF).I^^I^ 


We examine energies and double occupancies, which 
are single numbers and can be obtained by essentially all 
methods, enabling straightforward comparison. We also 
consider properties related to the electron Green’s func¬ 
tion, which at this stage are only available from a few 
methods. The results obtained from different techniques 
enables us to identify regimes of phase space that are 
well understood, in the sense that several different meth¬ 
ods provide results that are converged and agree within 
(reasonable) errors, and regimes that are not well un¬ 
derstood, in the sense that there is as yet no agreement 
between different methods. We show excellent agreement 
and small uncertainties between numerically exact tech¬ 
niques at half filling (all coupling strengths), weak cou¬ 
pling (all carrier concentrations) and for carrier concen¬ 
trations far from half filling (most interaction strengths). 
For carrier concentrations near to half filling and for in¬ 
termediate interaction strength, results can be obtained, 
but the resulting uncertainties are much larger in gen¬ 
eral, and more difficult to eliminate. We surmise that at 
least part of this uncertainty has a physical origin related 
to the presence of several competing phases, leading to 
sensitive dependence on parameters. 

The rest of this paper is organized as follows. In sec¬ 
tion]^ we define the model, delineate parameter regimes 
and define and discuss the observables of interest in this 
paper. In section [Hi] we summarize the methods, giving 
brief descriptions and focusing on issues most relevant 
to this paper while referring the reader to the literature 
for detailed descriptions. Section IV A demonstrates the 
importance of the extrapolation of results to the ther¬ 
modynamic limit and discusses the issues involved in the 
extrapolations. In sections |Vj |VI| and |VH| we present 
static observables in the strong coupling, intermediate 
coupling, and weak coupling regimes respectively. Sec¬ 
tion |VlTT] presents momentum and frequency dependence 
at finite T. A conclusion summarizes the work, out¬ 
lines the important areas where our present-day knowl¬ 
edge is inadequate, and indicates directions for future 
research. Supplementary material presents the thermo¬ 
dynamic limit values obtained here. A database of results 
is made available electronically.EsBsl 


II. THE HUBBARD MODEL 

The Hubbard model is defined by the Hamiltonian 
H = -Y^ Uj + H.c)j + G ^ (1) 

i 

where c\^ {ci^) creates (annihilates) an electron with spin 
a =t)i on site *, riia = is the number operator, 

and tij denotes the hopping term. In this work, we will 
restrict the discussion to the repulsive Hubbard model 
[U > 0) defined on a two dimensional square lattice; we 
further assume that the hopping contains only nearest 
neighbor {Uj = t) and second nearest neighbor {tij = t') 
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Parameters 

Parameters studied 

YJl 

-0.2 

0 

0.2 



U/t 

2 

4 

6 

8 

12 

n 

1.0 

0.875 

0.8 

0.6 

0.3 

T/t 

0 

1/8 

1/4 

1/2 



TABLE I. Parameters studied for the Hubbard model, t de¬ 
notes the nearest-neighbor hopping, t' the next nearest neigh¬ 
bor hopping, U the interaction strength, n the density, and T 
the temperature. 

terms, t is used to set the scale of all energies presented 
in this work. 

We consider interaction strengths ranging from U/t = 
2 to U/t = 12 and focus on temperatures where high- 
temperature expansion methods fail.^ We give repre¬ 
sentative results for the ground state and temperatures 
T/t = 0.125, 0.25, and 0.5. Table |l] contains a complete 
list of the parameters studied. 

At zero temperature, we systematically compute the 
energy per site and the double occupancy, and we present 
also some data on the nature of the order and the order 
parameter, where an ordered phase is found. At non-zero 
temperature, we also present dynamical information, in 
particular the Matsubara self-energy. We focus on val¬ 
ues at a given density, rather than at a given chemical 
potential. This implies that methods based on a grand- 
canonical formulation need to adjust the chemical poten¬ 
tial to find the right density, leading to additional uncer¬ 
tainty in computed quantities. 


III. METHODS 
A. Overview 

Many numerical methods provide solutions to the Hub¬ 
bard Hamiltonian on a finite size lattice. In this work, 
we restrict attention to techniques which can access sys¬ 
tems large enough that an extrapolation to the infinite 
system can be performed, and our aim is to obtain re¬ 
sults for the thermodynamic limit. Also important to our 
analysis is the assessment of uncertainties, either by pro¬ 
viding an unbiased error bar or an error bar that contains 
all errors except for a systematic contribution which may 
be assessed by comparison to other methods or reference 
systems. 

We consider three broad classes of methods: ground 
state wave function methods, embedding methods, and 
Green’s function methods. The distinction between these 
methods is not sharp; several of the algorithms fit into 
multiple categories, but the categorization provides a use¬ 
ful way to organize a discussion of the different kinds of 
uncertainties. 

Wave function methods construct an approximation to 
the ground state wave function of a system. Expectation 
values of observables (energies and correlation functions) 
are then evaluated by applying operators to this ground 


state wave function. The issues are the accuracy of the 
wave function for a given system size, and the accuracy of 
the extrapolation to the thermodynamic limit. AFQMC 
(with and without constrained path), UCCSD, FN, MR- 
PHF, and DMRG are wave function methods. 

Embedding methods approximate properties of the full 
system (for example the self energy or density matrix) by 
the solution of a finite cluster self-consistently embedded 
in an appropriately designed infinite lattice. The full so¬ 
lution of the original problem is recovered as the cluster 
size is taken to infinity. Errors in embedding methods 
arise from three sources: the solution of the cluster prob¬ 
lem, the convergence of the self-consistency loop that per¬ 
forms the embedding, and finite cluster size, with maxi¬ 
mum cluster sizes depending on the method and ranging 
from 16 to 100. The DMET, DCA, DF are embedding 
methods. 

Green’s function methods are defined here as meth¬ 
ods based on stochastic evaluations of many-body per¬ 
turbation series. They provide many-body self-energies, 
and Green’s functions, as functions 

of momenta, k, and fermionic matsubara frequency, itOrn 
from which other observables (energies and densities) can 
be calculated. Techniques which produce real-freauency 
(rather than matsubara frequency) informatioiP™^ are 
either restricted to small systems, molecules or impurity 
models or work best at weak coupling. DiagMC is a 
Green’s function method formulated directly in the in¬ 
finite lattice; the main issues for this method are the 
accuracy of the stochastic approximation to the full dia¬ 
grammatic expansion and the systematic truncation and 
extrapolation of the series. The DF and DGA techniques 
use Green’s function techniques to evaluate the impu¬ 
rity problem and the expansion around it; they therefore 
are subject both to embedding uncertainties and to the 
uncertainties arising from the evaluation of the diagram¬ 
matic expansion. 


B. Auxiliary-field Quantum Monte Carlo(AFQMC) 

Auxiliary-field Quantum Monte Garlo (AFQMG) is a 
ground-state wave-function method based on the idea 
that in the limit /3 —>■ oo the operator e~^^ applied to 
an initial wave function projects out the ground 

state of the Hamiltonian H. The projection is formulated 
as an imaginary-time path integral that is stochastically 
evaluated with the help of auxiliary fields introduced by 
a Hubbard-Stratonovich transformation. The method is 
applied to finite size lattices and an extrapolation to the 
infinite lattice case is required. If the calculation is con¬ 
verged, the exact ground state energy and wave function 
for the lattice are obtained. The issues are the conver¬ 
gence of the stochastic evaluation of the projector and, 
when particle-hole symmetry is broken, the presence of 
a sign problem. The sign problem is managed using a 
constrained path approximation, which introduces a po¬ 
tential systematic error that must be quantified by com- 
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parison to other methods. For an introduction to the 
basics of AFQMC methods, see, e.g., Ref. 

In this manuscript we present results obtained from 
two ground-state auxiliary-field quantum Monte Carlo 
(AFQMC) methods. The firslP^ is based on the ground- 
state path integral form of AFQMC but intro¬ 
duces several new algorithmic ingredients inclu ding an 
acceleration techniqu^^ (with force bia^^^^ in the 
Metropolis sampling and control of Monte Carlo variance 
divergence.l^This approach is applied to wstems at half 
filling with f = 0, where the sign problerrP^ is absent be¬ 
cause of particle-hole symmetry.l^ The algorithmic im¬ 
provements allow us to reach longer imaginary time in 
the calculations, achieve a higher acceptance ratio, and 
greatly reduce the Monte Carlo variance.l^ The second 
approach we employ, to treat cases where the sign prob¬ 
lem does occur, is re ferred to as the constrained path 
Monte Carlo method.^^^^^ This approach controls the 
sign problem with a constraint (implemented via a choice 
of trial wave functiorP^^ on the paths in auxiliary- 
field space, which allows stable calculations for arbitrarily 
long imaginary time and system size. 

Both methods obtain the ground state of the Hamil¬ 
tonian for a supercel l of si ze L x L under twist-averaged 
boundary conditions The ground state is obtained 

by use of Monte Carlo methods to estimate = 

g-/3ff.0(/3=o))^ The total projection length is typically 
/I = 64 in the ground state projection method, al¬ 
though test calculations with imaginary-time lengths sev¬ 
eral times larger were performed. The convergence error 
from finite values of /3 is negligible. In the constrained 
path method, the runs are open-ended and tend to cor¬ 
respond to much larger values of /3. In both the ground 
state projection and constrained path methods the sta¬ 
tistical error from the Monte Carlo calculation can be 
reliably estimated (one-standard deviation error bars re¬ 
ported) . The systematic error from the constrained path 
method is not variational, but depends on |4'7’) in the 
sense that it vanishes if I'I't’) is exact. Its magnitude will 
be quantified below by comparison to other techniques. 

A Trotter decomposition is used in the imaginary time 
evolution. The Trotter error from the finite time step. 
At = P/n, must be extrapolated to zero. This extrapo¬ 
lation can be controlled and does not make a significant 
contribution to the error budget. Most calculations re¬ 
ported here use a time-step fixed at Ar = 0.01 in units 
of t. 

Results obtained for finite L are averaged over twist- 
angle, 0, to remove one-body finite size effects. For small 
systems a large nu mber of 0 values are needed (~ 200 
for 4 X 4 or 6 X 6)1111111111 but for larger systems, far 
fewer 0 are required to reach the same level of accu¬ 
racy (for a L X L system with L = 20, averaging over 5 
twist angles is sufficient). These results are then extrap¬ 
olated to L —>■ oo. The extrapolation requires care be¬ 
cause the ground state depends on the system geometry. 
For n = 0.875, we used rectangular supercells (mostly 
8 X 32, checked with sizes 8 x 16, 16 x 16 and 8 x 48 for 


consistency) to accommodate spin- and charge-density 
wave orders. The extrapolation also requires careful at¬ 
tention to the functional form of the leading finite volume 
correction.lSSHlll Qur final results at the thermodynamic 
limit include all statistical errors, and a conservative esti¬ 
mate of the uncertainty resulting from the extrapolation 
of L —> oo in order to remove the two-body finite-size ef¬ 
fects (The fit includes 1/L^ and 1/L^ term^^IHlSl for the 
energies, and 1/L for magnetization m^). To provide a 
concrete example, at n = 1, t' = 0, T = 0, ?7/t = 4, the 
energy per site for a 20 x 20 system after 0-averaging 
is E/t = —0.86038 ± 3 x 10“®. After a weighted least 
square fitting with L = 4, 6, • • • ,18, 20, the final result in 
the thermodynamic limit is E/t = —0.8603 ± 2 x 10“'^. 

For n ^ 1 or 7 ^ 0, a sign problem appears. The 
sign problem makes it impossible to converge the ground 
state projection method for the system sizes and propa¬ 
gation lengths /3 needed and an alternative method, the 
constrained path (CP) approximation is used. The re¬ 
sults reported in this paper follow Ref. HT] using a trial 
wave function jd^r) to apply a boundary/gauge condi¬ 
tion on the paths that are included in the path integral 
in auxiliary field space. All results reported here used 
single-determinant IT^-) with no release. In these cal¬ 
culations, I'I't’) is taken to be a mean-field solution for 
the Hubbard model for given U, L and 0 with a U value 
UeB = min{?7, At}. The order parameter in the mean-field 
solution is chosen to be orthogonal to the spin quantiza¬ 
tion axis. This choice is found to hel p pre serve symme¬ 
try in I'kj’), improving the CP results.^^SlSIl xbe accuracy 
of the constr ained path approach has been extensively 
benchmarked.l22ISiEIEll Pave carried out additional 
comparisons with exact diagonalization on 4 x 4 systems. 
At U = At, the relative error on the energy, averaged over 
60 randomly chosen 0 values, is -1-0.018% for n = 0.25 
and -0.15% for n = 0.625. At U/t = 8 and n = 0.875 
the relative error is -0.51% averaged over 20 random 0 
values. We have also verified these estimates in a few 
systems of larger L, using mult i-dete rminant trial wave 
functions and constraint release. 


C. Fixed node Diffusion Monte Carlo (FN) with 
nodes from Variational Monte Carlo (VMC) 

The variational Monte Carlo method constructs a trial 
wave function that approximates the exact ground state 
of a correlated Hamiltonian.l22H22l xpe wave function de¬ 
pends on parameters that are optimized by minimizing 
the expectation value of the Hamiltonian, which requires 
a Monte Carlo sampling whenever the trial state is cor¬ 
related (i.e., it is not a simple Slater determinant). We 
remark that the variational Monte Carlo energy gives an 
upper bound to the exact value and that, with this ap¬ 
proach, it is possible to access quite large clusters, with 
all relevant spatial symmetries (translations, rotations, 
and reflections) preserved. However, it is difficult to 
quantify the systematic errors introduced by the choice of 
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the trial state. Moreover, while spatial correlations may 
be correctly captured, dynamical properties are missed. 

We generated simple variational wave functions by ap¬ 
plying a density Jastrow factor on top of uncorrelated 
states that are built from local (mean-field) Hamiltoni¬ 
ans, containing only few parameters, where the physical 
properties are reflected in a transpa rent way as differ¬ 
ent terms inside the variational state.^^^*^ At finite dop¬ 
ing, the uncorrelated states have been obtained from the 
BCS Hamiltonian, including superconducting pairing on 
top of electron hopping; in addition, collinear antiferro¬ 
magnetism with Neel order parallel to the 2 spin quan¬ 
tization axis is also included. At half filling, where the 
system exhibits long-range magnetic order, the uncorre¬ 
lated state contains only magnetism in the x — y plane; 
in this case, an additional spin Jastrow factor involv¬ 
ing the ^-component of the spin operator is also taken. 
This term couples spins along a direction orthogonal to 
the magnetic ordering plane, reproduciM the spin-wave 
fluctuations above the mean-field state.^^Sl All these vari¬ 
ational wave functions with Jastrow factors generalize 
the so-called Gutzwiller state, ^ allowing a descripti on of 
metals, superconductors, and also of Mott insulators.^^^^^ 
Nevertheless, they do not give an accurate approximation 
to the exact ground state in two spatial dimensions, es¬ 
pecially close to half filling. We obtain a substantial im¬ 
provement by including the backflow correlations inside 
the uncorrelated state. On the lattice, this corresponds 
to a redefinition of the single-particle orbita ls and it is 
particularly efficient at strong couplingj ^^ l ^^ l 

To determine the variational parameters (i.e., the ones 
related to the Jastrow factors, the ones included in the 
mean-field Hamiltonian, and also those of the backflow 
correlations), we minimize the expectation value of the 
Hamiltonian. This minimization is performed by con¬ 
structing a Markov chain using the Metropolis algorithm, 
where walkers are defined by many-body configurations 
having electrons on lattice sites with given spin along 
the z axis. After performing this optimization, a further 
improvement can be obtained by applying the Green’s 
function Monte Carlo projection technique^ to the opti¬ 
mal trial state within a fixed-node approximation.!^ This 
procedure allows accurate calculations of the energy and 
diagonal correlation functions, such double occupancies 
or density-density correlation functions. The Ansatz on 
the nodal structure given by the variational wave func¬ 
tion induces a systematic error, which cannot be deter¬ 
mined a priori but can be estimated from the change in 
energy as the trial wave function is improved. We point 
out that there is a difference between continuum and lat¬ 
tice fixed-node approaches. In the continuum, only the 
signs of the trial function are important: if the nodes 
are correctly placed, the exact energy is obtained. By 
contrast, on the lattice, both the signs and the relative 
magnitudes of the trial function in configurations that 
are connected by a sim flip must be correct in order to 
have the exact energy.1221 The error bars for finite systems 
are given as the statistical errors of the Green’s function 


Monte Garlo technique and do not include any estimates 
of the systematic errors coming from the fixed-node ap¬ 
proximation. 

The finite size results are then extrapolated to the in¬ 
finite system size by using a scaling that depends on the 
carrier concentration. At half filling, we use the 
scaling (where N is the system size) that is appropri¬ 
ate for two-dimensional ordered antiferromagnets.!s3^ In 
this regime, the error bar for the infinite system size is 
given by a fitting error of the linear regression. At finite 
dopings, the size scaling may suffer from shell effects: a 
smooth behavior can be obtained only when a sequence 
of closed-shell configurations are taken (i.e., electron fill¬ 
ings for which the non-interacting case corresponds to a 
unique ground state). In the generic case, size effects may 
be dominated by the ones present at t/ = 0. This is the 
case for large dopings (e.g., n = 0.8) and all interactions, 
and intermediate dopings (e.g., n = 0.875) and small 
interactions {U < 4). Here, for every available size, the 
ratio between the energy at finite U and the one at ?7 = 0 
is roughly constant and the thermodynamic limit can be 
assessed by fitting this ratio, namely the infinite-size en¬ 
ergy is obtained by multiplying the aforementioned ratio 
by the thermodynamic value at t/ = 0. The extrapolated 
value is assumed to be normally distributed with an er¬ 
ror bar taken as the difference between the estimated 
thermodynamic limit and the largest available size. For 
intermediate dopings (n = 0.875), the size scaling starts 
to deviate from the t/ = 0 case around t//t = 4, and we 
decided to take the point at the largest size as the infi¬ 
nite size limit, with an uncertainty of twice the difference 
to the next lowest system size. We remark that, in this 
case, a linear regression with the 1/A^ scaling gives an 
estimate of the thermodynamic limit that is compatible 
(within one error bar) with the point at the largest size. 


D. Multi-reference Projected Hartree—Fock 
(MRPHF) 

The mult i-reference projected Hartree-Fock 

methocPSllillll 

is a ground state wave function ap¬ 
proach based on a trial wave function jT) characterized 
by the quantum numbers 0, K that is constructed out of 
a set of broken-symmetry Hartree-Fock wave functions 
via projection operators. The idea is that a 
broken-symmetry determinant includes the dominant 
correlation physics while the projection restores the 
physical symmetries. The wave function takes the form 

n 

\'^K)=T.T.fK^PKK'm, (2) 

i=l K' 

and the parameters are determined by minimizing 
the energy. The projector P^k' restores the symmetries 
(characterized by the quantum numbers 0, K) in j^*) and 
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can be formally written a^^^IZSl 

Prk' = ^ ^ R{m) (3) 

ra 

in terms of the rotation operators R{m) and the irre¬ 
ducible representation matrices r®(TO) associated with 
the elements m of the symmetry group of the problem. 
Here, h is the dimension of the irreducible representation, 
while L is the volume of the group. The character of the 
broken symmetry determinant is optimized in the pres¬ 
ence of the projection operator (i.e., a variation-after- 
projection approach), which results in broken symmetry 
determinants with well defined defects.!^ 

Our expansion employs Slater determinants that break 
the space group and spin {S'^) symmetries of the lattice, 
but preserve Sz symmetry. All the broken symmetries are 
restored using the appropriate projection operators. The 
series of i determinants in Eq. is constructed through 
a chain of variational c alculatio ns, using the FED (FEw- 
Determinant) approachIn this procedure, after a 
wave function with n — 1 intrinsic determinants is avail¬ 
able, a wave function with n determinants is variationally 
optimized by adjusting the Thouless coefficients deter¬ 
mining the last added determinant. The full set of linear 
coefficients is re-adjusted. The MR-PHF approach 
becomes exact as the number of determinants retained 
tends to the size of the Hilbert space, i.e., as Eq. Q be¬ 
comes a coherent state representation of the exact ground 
state wave function. 

If the number of determinants is fixed at a finite, not 
too large, value, these calculations can be performed for 
large lattices with polynomial cost {0{N)^ or so, where 
N is the number of sites in the lattice). However, the 
number of determinants required to obtain results with 
a given accuracy increases exponentially with increasing 
system size. In this work, implementation aspects have 
compelled us to keep the number of determinants roughly 
constant, so that the quality of the solution decreases as 
system size is increased, and, consequently, the energy in¬ 
creases, precluding a thermodynamic limit extrapolation. 
We have used expansions with 4, 24, and 32 determinants 
for half-filled 4 x 4, 6 x 6, and 8x8 lattices, respectively. 
In the lightly-doped regime, larger expansions have been 
used: 48 determinants in a 10 x 4 ((n) = 0.8) lattice and 
80 determinants in a 16 x 4 ((n) = 0.875) one. A linear 
extrapolation in the reciprocal of the number of deter¬ 
minants, i.e. in 1/n has been performed to the infinite 
configuration limit for the ground state energies. 

The calculations presented could in principle be im¬ 
proved in a number of ways: additional symmetries 
could be broken and restored (such as Sz or particle 
number) in the reference configurations; more config¬ 
urations could be included; and/or a full optimization 
of all determinants could be performed, in the spirit of 
the resonating Hartree-Fock approach.^S The accuracy 
of any one result can therefore be increased, as shown W 
Rodrfguez-Guzman et al^ or by Mizusaki and ImadeP^ 


in the closely related path integral renormalizati on g roup 
(PIRG) approach. Recently, Tahara and ImadaP^ have 
combined the symmetry-projected determinant expan¬ 
sion with short- and long-range Jastrow factors within a 
variational Monte Carlo framework, which may be used 
to further increase the accuracy. These techniques, and 
others, have been used to explore, for example, spin and 
char ge s tripe phases, which we do not explore in this 

workPHU] 


E. Unrestricted Coupled Cluster - Singles and 
Doubles (UCCSD) 

Coupled cluster (CC) theorjBsHUl is a ground-state 
wave function technique. It is widely used in quantum 
chemistry and often considered the best source of pre¬ 
cise data for molecules that are neither too large nor 
too strongly correlated. Its application to the Hubbard 
model has been more limited, where it has been used in 
two different formulations. In the first form, which di¬ 
rectly exploits the translational invariance of the lattice 
to work in the thermodynamic limit, the theory is formu¬ 
lated in the site-basis, starting from an infinite Neel or¬ 
dered reference, from which clusters of excit ation s which 
change occupancy and flip spins are createcP^MI jn the 
second form, the theory starts from a single determinant 
reference state on a finite lattice, from which clusters of 
particle hole excitations are createcP^. This is similar to 
how the theory is used in quantum chemistry, and is the 
formulation discussed further here. As this second form 
does not work in the thermodynamic limit, the energies 
must be extrapolated. 

The GG wave function is written as |'I') = exp(T) |0) 
where |0) is a single determinant reference state, and 
T = J2n'^n, where Tn = is file cluster op¬ 

erator. The operator creates an excited determi¬ 
nant |$q^) which contains n particle-hole pairs relative to 
the reference state. In its standard and simplest version 
(used here), the energy and coefficients tq^ are obtained 
by solving the Schrodinger equation projectively: 

E = (0|e-^i7e^|0), (4a) 

0 = (4>,^|e“^i7e^|0) Wqn- (4b) 

CG theory thus diagonalizes a similarity-transformed 
Hamiltonian H = exp(—T) H exp(T) in a subspace of 
states defined by |0) and 1$,^). Note that because T 
is a pure excitation operator, the commutator expansion 
used to evaluate H truncates after four commutators. 

If the sum over n in defining the cluster operator is 
carried out to all orders, the exact ground state wave 
function is reproduced. In practice the operator T is 
typically restricted to terms involving a small number 
of particle-hole pair excitations above the reference state 
(n < Umax) and l^*) = exp(T) |0) is projected onto the 
space of states with up to Umax particle-hole pairs; the 
accuracy then depends upon nmax and the choice of ref¬ 
erence state |0). In a lattice model such as the Hubbard 
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model with N sites, the computational cost is, roughly 
speaking, proportional to 

The calculations reported in this manuscript primarily 
limit T to n < 2, i.e. to the creation of only singly- and 
doubly-excited determinants, giving what we r efer to as 
CC with single and double excitations (CCSD)!^^*^ For 
select examples, we have included corrections for triple 
and occasionally quadruple excitation effects. The accu¬ 
racy of CC theory, and the need for higher excitations, 
depends on how well the reference |0) captures the quali¬ 
tative physics. When the reference is accurate, then sin¬ 
gles and doubles excitations may be sufficient; however, 
when the reference determinant bears little resemblance 
to the exact wave function (as may happen in strongly 
correlated systems) a much higher degree of excitation is 
required to recover the correct physics. For this reason, 
the calculations reported here use a symmetry-broken 
unrestricted Hartree-Fock (UFfF) reference determinant, 
because UHF can provide a better mean-field descrip¬ 
tion, particularly near half filling where antiferromag¬ 
netic correlations dominate; this defines the unrestricted 
CCSD (UCCSD) method used here. Note, however, that 
particularly for doped systems with large U there are 
a plethora of nearly degenerate UHF states, and find¬ 
ing the best reference for UCCSD is not straightforward. 
We have prepared the UHF solution following the pre¬ 
scription of Ref. |53 In principle the deficiencies of the 
reference determinant can be corrected in what is known 
as Brueckner ccP where the reference determinant is 
adjusted to eliminate single excitation effects. These cal¬ 
culations are more computationally demanding and we 
have not pursued them here. 

An important virtue of the exponential parametriza- 
tion of the wave function is that the CC energy has a 
non-trivial thermodynamic limit even for restricted ex¬ 
citations rimax^^^ Thus, as lattice size increases, the 
energy per site approaches the thermodynamic limit for 
the given Umaxi S'lid the exact thermodynamic limit as 
nmax is increased. For smaller values of U where conver¬ 
gence to the thermodynamic limit is slower, we have con¬ 
verged second-order perturbation theory out to the ther¬ 
modynamic limit and added a correction for the differ¬ 
ence between CC theory and perturbation theory which, 
for small C/, converges quickly with respect to system 
size. Double occupancies have been computed by nu¬ 
merical differentiation of the CC energy with respect to 
U. We also provide UCCSDTQ estimates of the ground 
state energy (labelled as UCCSDTQ*) for n=l systems, 
where the triples (T) correction is obtained as UCCSDT- 
UCCSD energies for a 6 x 6 system, and the quadruples 
(Q) correction is obtained from UCCSDTQ-UCCSDT 
energies for a 4 x 4 lattice. For n=0.8 (n=0.875), our 
UCCSDT estimates (labeled as UCCSDT*) are obtained 
from 10 X 4 (16 x 4) lattices. 

The UCCSD calculations reported here can be com¬ 
pleted in a few hours using standard quantum chemistry 
packages.Even at half filling where there are not many 
Hartree-Fock solutions to be concerned with, we find 


large effects due to single excitations, which suggests that 
the coupled cluster calculations could be substantially 
improved by optimizing the identity of the reference de¬ 
terminant. Similarly, we generally see significant correc¬ 
tions due to triple and higher excitations; these can also 
be computed with standard packages,® but while opti¬ 
mizing the reference determinant and including higher 
excitation effects will increase the accuracy of the cou¬ 
pled cluster calculations, they also increase the cost. 


F. Density Matrix Renormalization Group 
(DMRG) 

The density matrix renormalization groujJ® is a vari¬ 
ational ground state wave function technique. It con¬ 
structs the ground state of a system by diagonalizing the 
Hamiltonian in a finite subspace spanned by an itera¬ 
tively constructed basis that is optimized via a Schmidt 
decomposition that minimizes the spatial extent of the 
quantum mechanical entanglement between basis states. 

DMRG is generally believed to be the optimal method 
for finding grounds states of one dimensional lattice mod¬ 
els. In the application to one dimensional systems two 
sources of error must in principle be controlled. Results 
for a fixed system length of L must be converged with 
respect to basis size m, and then the converged results 
must be extrapolated to L —>■ oo. However, for most of 
the one dimensional Hamiltonians of current interest the 
convergence is very rapid and in practice large enough m 
and L are accessible numerically, so extrapolation is not 
required. 

Application of DMRG to a finite size 2D system pro¬ 
ceeds by defining an effective one dimensional problem, 
to which the standard one-dimensional DMRG is applied. 
Two dimensional tensor network generalizations of the 
DMRG ideas have attracted tremendous recent interest, 
but these methods have not yet produced results for the 
two dimensional Hubba rd mo del that can be included in 

the present comparison.®^!] 

Most current implementations of DMRG require open 
boundary conditions. The canonical method for creat¬ 
ing an effective one dimensional system from a finite-size 
two dimensional one is to impose periodic boundary con¬ 
ditions in one direction and open boundary conditions 
in the other, thereby defining a cylinder of finite length 
and finite circumference. One then defines an effective 
one dimensional problem by indexing the sites along a 
one dimensional path that covers all of the sites on the 
cylinder and taking the matrix elements of the Hamil¬ 
tonian in this basis. The price is that two sites sepa¬ 
rated by a small distance along the cylinder axis in the 
physical system are separated by a distance of order the 
cylinder circumference in the effective one dimensional 
model. The effective one-dimensional problem thus has 
long ranged terms, which imply longer ranged entangle¬ 
ment and require that more states are kept in the optimal 
basis. The number of states needed grows exponentially 






with the circumference of the cylinder (width of origi¬ 
nal finite lattice) meaning that there is a sharp cut-off in 
the accessible system widths, typically around width 6 
in the Hubbard systems; however, systems of very large 
cylinder length can be studied. The extrapolation to the 
thermodynamic limit must thus be handled with care. 

The DMRG calculations reported here were perfo rmed 
with the standard DMRG finite system algorithnPSEIl 
Two types of cylinders were considered: one with the ax¬ 
ial and circumferential directions aligned parallel to the 
bonds of the square lattice and one rotated by 45 degrees. 
When one cuts a cylinder in two, the 45 degree rotated 
system has fewer sites on the boundary per unit length, 
and thus one expects a smaller growth of the entangle¬ 
ment with the length of the cut (governed by the area 
law). For the undoped antiferromagnetic system, the ro¬ 
tated system also is unfrustrated both for odd and even 
circumferences, reducing shell effects in the finite size re¬ 
sults. For half filled systems both types of orientation 
are considered and they show good agreement, and error 
bars are estimated to incorporate both the error bars on 
the data points for specific widths, and the differences be¬ 
tween the two orientations. However we note that for the 
half-filled syste ms, b etter results were obtained with the 
rotated system.l22Ely\tith doped systems, we see striping 
behavior at stronger coupling. A stripe is a line of holes 
which act as a domain wall in the antiferromagnetism 
on either side. These stripes have been seen in the t-J 
model with DMRG starting in the late 90’s.l^The stripes 
typically wrap around the cylinder, with a specific even 
number of holes in the ring-stripe. With doping, the op¬ 
timal number of holes can change. Striped configurations 
with the wrong number of holes in a stripe are typically 
metastable. For the ordinary orientation, we were able 
to sort out stripe fillings, finding the low energy states 
and avoiding metastability. For the 45 degree rotated 
lattices, the patterns seem more delicate, and we have 
not yet sorted out the lowest energy configurations. We 
thus present only the results for standard orientation for 
doped systems. 

Gonvergence issues pose more severe problems than in 
the standard one dimensional cases, both because of the 
intrinsic limits on system size discussed above and be¬ 
cause in some cases the presence of several metastable 
states can cause troubles for extrapolation, and can lead 
to the appearance of states that are not the ground state 
and may be important only for finite systems. In the sys¬ 
tems studied here, metastability was traced to the pres¬ 
ence of physically different “striped” states for different 
hole doping in the DMRG cylinder. 

We obtained converged results as follows. For each 
cylinder of a specific length L, we extrapolated the energy 
in the number of basis states, m. To make this extrapola¬ 
tion reliable m was slowly increased, but each m was used 
for two consecutive sweeps. The truncation error and the 
energy were measured on the second sweep for each m. 
Then a linear extrapolation of energy versus truncation 
error is used to obtain the ground state energy with error 


bars. The deterministic nature of DMRG can result in 
uncertainties due to fitting which do not appropriately 
represent the uncertainty in choice of extrapolation pro¬ 
cedure. We therefore assume a normally distributed error 
of 1/5 the difference between the last point and the ex¬ 
trapolated value which we justif y by c omparison to the 
accuracy of previous DMRG data.l2212il Metastability was 
signaled by lack of linearity in this extrapolation. When 
this was found we determined which state had the lowest 
energy. The system was then rerun with an initial state 
favoring the lower energy state producing results in the 
lower energy configuration such that the extrapolations 
are linear in the truncation error. 

For a fixed width we then extrapolate the energies lin¬ 
early in 1/L to get an energy per site for an infinite 
cylinder. Errors are estimated statistically using the er¬ 
ror bars on each point. To reduce the finite size effects 
from different widths, we employ a simple version of the 
phase averaging triclP^ by taking an average over pe¬ 
riodic and antiperiodic boundary conditions around the 
cylinder. The phase averaging eliminated an oscillation 
in the energy as a function of system width in the 45 
degree rotated systems, and in general nicely accelerated 
convergence. We then analyzed the results as a func¬ 
tion of cylinder width. For half filling we found that 
extrapolating the energy to the thermodynamic limit by 
1/width^ works well. This is the finite-size behavior in 
the Heisenberg model, where it is well understood.!^ 


G. Density Matrix Embedding Theory (DMET) 

The density matrix embedding theory (DMET)pSl^is 
a ground state embedding method formulated in terms 
of wave function entanglement. Given an impurity clus¬ 
ter of size Ac, DMET maps an L x L bulk many-body 
problem (L is chosen very large) to an impurity and bath 
many-body problem, yielding a problem with 2Nc sites 
in total. The mapping is constructed from the Schmidt 
decompositiorP^ of a bulk wave function j^/). The formu¬ 
lation is exact if I'k) is the exact bulk ground state, or in 
the limit of impurity cluster size Nc —> oo. Since the ex¬ 
act bulk state is not known a priori, an approximate bulk 
state is used for the impurity mapping. Recently, DMET 
has been applied to both ground-state and linear re¬ 
sponse spectral properties of the Hubbard modePsEHMl 
In this work we use a general BCS bulk state, the ground- 
state of the DMET lattice Hamiltonian given by the hop¬ 
ping part of the Hubbard Hamiltonian augmented with 
the DMET correlation potential u, that is applied in each 
cluster supercell of Nc sites in the bulk lattice. This bulk 
state is allowed to spontaneously break spin and parti¬ 
cle number symmetry through the self-consistency cycle 
that determines u (which contains both particle number 
conserving and non-conserving terms). In this cycle, the 
bulk state |4t) is updated from the interacting impurity 
and bath solution 1$), by minimizing the difference be¬ 
tween l^/) and 1$) (as measured by their (generalized) 
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one-particle density matrices) with respect to the poten¬ 
tial u. 

For the bulk lattice, we use L = 72. From cal¬ 
culations on larger L, we find that the finite lattice 
error associated with this choice is negligible, on the 
scale of the significant digits reported. The BCS bulk 
state is obtained by solving th e lattic e spin-unrestricted 
Bogoliubov-deGennes equatiorfi^Slini] -with the correla¬ 
tion potential u. The impurity and bath problem is 
solved in the BCS quasiparticle basis, with general one- 
body and two-body interactions that do not conserve 
particle number or locality. We use the density ma¬ 
trix renormalization group (DMRG)I^ as an impurity 
solver (adapt ed from the quantum chemistry DMRG 
code BLOCKliSMloil)^ with a maximum number of renor¬ 
malized states m=2000 (DMET self-consistency is per¬ 
formed up to m=1200). The DMET lattice and impurity 
Hamiltonians are augmented with a chemical potential 
/r, to ensure that the relative error in particle number is 
less than 0.05%. We used impurity clusters of dimensions 
2x2, 4x2, 4x4 and 8 x 2 at each point in the phase 
diagram. The energies and observables were then extrap¬ 
olated to the thermodynamic limit using a linear relation- 
ship with Nc ' , as appropriate to a non-translationally 
invariant cluster embedding theory. 

The total DMET uncertainty is estimated by com¬ 
bining the errors from three sources (i) convergence 
of DMET self-consistency, (ii) solution of the impurity 
many-body problem using DMRG, (iii) extrapolation to 
the limit of infinite impurity size. We estimate the 
self-consistency error (i) using the difference of the last 
two DMET self-consistent iterations. The average self- 
consistency error is below 5 x in the energy for all 

cluster sizes. The impurity solver error (ii) is from us¬ 
ing a finite number of renormalized states m in DMRG. 
This error is only non-zero in clusters larger than 2x2. 
The energy and observables are extrapolated to to = oo 
using the standard line ar relatio n between energy and 
DMRG truncation error.l2 ^°^ l ^°^ l For 4x4 impurity clus¬ 
ters, the truncation error is large enough to contribute 
also to the converged DMET self-consistent correlation 
potential u{m). To take this into account, we extrap¬ 
olate using (a) self-consistent DMET results converged 
at different to, and (b) non-self-consistent DMET results 
using different to in the DMRG impurity solver at a fixed 
correlation potential u(TOniax)) where TOmax is the max¬ 
imum TO used in the DMET self-consistency. The dif¬ 
ference between the two extrapolations is then added to 
the total DMET uncertainty. For the cluster sizes in this 
study the errors due to (i) and (ii) are small and eas¬ 
ily controlled. Therefore, the finite size impurity error 
is the main source of uncertainty at almost all points in 
the phase diagram. It is estimated as the standard devi¬ 
ation of the finite size extrapolation. The quality of the 
approximate DMET impurity mapping depends on the 
approximate lattice wave function j'k) and decreases as 
the coupling strength U increases, especially for carrier 
concentrations near half filling. This slows down the clus¬ 


ter size convergence of the results for large U, increasing 
the uncertainty. In the strong coupling, weakly doped 
region, we find competing homogeneous and inhomoge¬ 
neous orders that become very sensitive to the cluster size 
and shapes (similar to the stripes observed in DMRG). 
It is difficult to reliably extrapolate these results to the 
thermodynamic limit. As a result, the total DMET un¬ 
certainties range from about at half filling and low 

densities to a maximum of about I0“^f in the strong cou¬ 
pling, underdoped region. A detailed description of the 
methodology and extrapolation procedures for the calcu¬ 
lations is contained in Ref. 11061 


H. Dynamical Cluster Approximation (DCA) 


The dynamical cluster approximation (DCA)I ^F07 | i08 | 
is an embedding technique in which an approximation to 
the electron self energy is obtained from the solution of a 
quantum impurity model consisting of some number Nc 
of interacting sites coupled to an infinite bath of non¬ 
interacting electrons. In applications to the Hubbard 
model, the interactions in the impurity model are the 
interactions of the original problem while the one-body 
terms are determined from a self-consistency condition 
relating the Greens functions of the impurity model to 
those of the lattice. The DCA is a particular general¬ 
ization to Nc > 1 of the ‘single site’ dynamical mean 
field method^^E^ Other generalizations have also been 
introduced, F^^ I ^^^ but results from these methods are not 
considered here. The single-site method was motivated 
by the observation that in an appropriately defined infi¬ 
nite coordination number limit an exact solution of the 
Hubbard model can be found,^and can be reca st in terms 
of the solution of a single-site impurity model.^^^ It was 
later understood that generalization to Nc > 1 impurity 
model representations allows a treatment of the finite 
coordination number model that becomes exact as the 
number of impurity sites Nc —>■ oo. 

The DCA formulation partitions the Brillouin zone 
into Nc equal area tiles and approximates the self en¬ 
ergy E as a piecewise constant function of momentum 
taking a different value in each tile: 

E(k,cc)= ^ </>a(k)E,(a;) (5) 

a=l...N^ 


with ^a(k) = 1 for k S a and (f>a0^) = 0 for k ^ a. 
The tiles, a, map directly on to impurity model single¬ 
particle levels, the impurity model Greens function and 
self energy are diagonal in the impurity model o-basis 
and the DCA self-consistency equation 


Gr^(cc) 


f cPk I 

a (27r)2 UJ - Ek- Ea(w) 


( 6 ) 


determines the remaining parameters of the impurity 
model. The self consistency loop is solved by iteration; 
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an initial guess for the impurity model parameters pro¬ 
duces a set of Ea which are then used to update the im¬ 
purity model parameters. The loop typically converges 
well and errors associated with the self consistency are 
smaller than errors in the solution of the impurity model. 

We obtain results for different W in the paramagnetic 
phase and extrapolate to the thermodynamic limit by 
exploiting the knowrP^ Nc scaling for momentum av¬ 
eraged qua ntities i n d-dimensions and systematically in¬ 
creasing 

To solve the W site impurity problem we u se a CT- 
AUX algorithnJS^E^ with submatrix updates Our 
codes are based on the ALPS libraries^^'^^ Selected 
points have been compared to a C T-IN frf^^ implementa¬ 
tion based on the TRIQS libraries .^1^ 

In this work we provide extrapolated DCA results from 
clusters of sizes W=16, 20, 32, 34, 50, 64, 72 and 98, 
depending on temperature and densities, in order to give 
a reliable estimate of the properties of the 2D Hubbard 
model in the thermodynamic limit. 

The CT-AUX method is a type of diagrammatic Monte 
Carlo. In the T > 0, impurity model context the dia¬ 
grammatic series is absolutely convergent and the issues 
discussed below for Diag-MC in the infinite-lattice con¬ 
text are not important. However, the CT-AUX method 
is restricted at low T by the existence of a Monte Carlo 
sign problem in the auxiliary held solver. The sign prob¬ 
lem worsens rapidly as U is increased, as T is decreased, 
or as Nc is increased and is particularly evident in the 
density range n = [0.8,1.0). Further, as temperature is 
decreased, the length scale of correlations in the system 
increases, resulting in larger hnite size effects. We take 
a conservative approach to determining the uncertainty 
in the extrapolation procedure. We include both a sta¬ 
tistical uncertainty in the extrapolation in l/W as well 
as an additional uncertainty which we take as half the 
difference of the extrapolated value from the largest Nc 
value explored. This gives a robust representation of an 
extrapolation error which is larger when hnite size effects 
are large, but that also vanishes as W increases. As such, 
the error bars for extrapolation of our DCA results to the 
thermodynamic limit contain both stochastic and hnite 
size uncertainties, and values for hnite system sizes with 
stochastic error bars are provided. 


I. Dual Fermion (DF) Ladder Approximation 

The dual fermion approacI l^^b 22 |i 2 ^ jg diagrammatic 
extension of the single site dynamical mean held theory 
(DMFT). The DF technique is motivated by the idea that 
non-local corrections to DMFT can be captured by a per¬ 
turbative expansion around a solution of the dynamical 
mean held equations. In formal terms the expansion re¬ 
quires summing a series of diagrams for two and higher 
particle correlations, with vertices dehned in terms of 
the fully interacting but reducible vertices of the impu¬ 
rity model. In this regard, the DF technique is similar in 


spirit to methods such as the dynamical vertex approx¬ 
imation (UFA) and dynamical mean-held extensions of 
fRG (DMF^RG) which approxim ate inte ractions on the 
level of 2-particle vertex functions.l^^lHilllin practical im¬ 
plementations to date the dual fermion expansion is trun¬ 
cated at the two particle level (higher than two-particle 
interactions are omitted) and the series of two particle 
corrections is approximated by a few low-order terms or 
a ladder resummation. One of its strengths lies in the 
ability to describe phase transitions of the syste m by em - 
ploying resummations of the relevant diagrams .112111211 
The DF results presen ted here are obtained using the 
open-source opendf cod d^^^ l, starting from a single-site 
{Nc = 1) dynamical mean held solution ob tained w ith 
a continuous-time auxiliary-held (CT-AUX)ii51hA0 im¬ 
purity solver. The method is limited by the accuracy 
to which reducible impurity vertex functions can be ob¬ 
tained, which is a polynomial (cubic) complexity prob¬ 
lem. Within the approximation of neglecting higher or¬ 
der vertices and only considering ladder diagrams, a sys¬ 
tematic estimation of deviation from the true interacting 
system is not possible, and we omit error bars altogether. 


J. Diagrammatic Monte Carlo (DiagMC) 


The diagrammatic Monte Carlo method (DiagMC) be¬ 
gins from the observation that within standard many- 
body perturbation theory, any quantity Q that depends 
on some set of arguments y (which may include both con¬ 
tinuous components such as frequency and momentum 
and discrete components such as spin) may be expressed 
as an inhnite series of terms, each of which involves multi¬ 
dimensional integrals and sums over many variables: 


j -j dxi... dXaD{a,^,xi,... ,Xa;y) , 


a=0 5 


(7) 

Here the D are known functions defined by the Feyn¬ 
man rules of diagrammatic perturbation theory. The 
series order a controls the number of internal integra¬ 
tion/summation variables {xi,..., Xq,}, and ^ labels dif¬ 
ferent terms of the same order in the series. The quantity 
Q may be the electron Greens function G, the self-energy, 
E, the screened interaction, W, the polarization opera¬ 
tor, H, the pair-propagator (for contact interaction), F, 
and its self-energy, Er, or the current-current or other 
correlation functions. Basic thermodynamic properties 
(particle, kinetic, and potential energy densities) in the 
Grand Canonical ensemble are readily obtained from G 
and E, see Ref. 11301 

The most widely-used formulation of perturbation the¬ 
ory is in terms of Feynman diagrams. In this case stan¬ 
dard rules relate the graphical representation of a given 
term in the series to the corresponding mathematical ex¬ 
pressions for the D, which are typically given (up to a 
sign or phase factor) by a product of functions associ¬ 
ated with graph lines, D = Hunes-^iine(xiine)- In pertur- 
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bation theory for particles interacting via pairwise forces, 
these lines are associated with the interaction potential, 
U, and bare single-particle propagators If we de¬ 

note the collection of all external and internal variables 
that allow for a complete characterization of the diagram 
(diagram topology and internal and external variables) as 
1 / = (a,^,xi,... ,Xa;y), then Eq. 0 can be viewed as 
weighted average over the configuration space {v}-. i.e., 
Q = where the modulus of D^, 

defines the configuration ’weight’, and ipu = &vgD^. The 
basic idea of the diagrammatic Monte Carlo (DiagMC) 
method is to use stochastic sampling of the configuration 
space to compute Q by treating oc \D,y\ as the probabil¬ 
ity density for the contribution of point v to the sum. In 
DiagMC the diagram order, its topology, and all internal 
and external variables are treated on equal footing and 
each diagram represents a point, not an integral, in {v}. 
The MC process of generating diagrams with probabili¬ 
ties proportional to their weight is based on the conven¬ 
tional Markov-chain updating schem^^lHlSl implemented 
directly in the space of continuous variables. 

Since connected Feynman diagrams are formulated di¬ 
rectly in the asymptotic limit, there is no infinite sys¬ 
tem size limit to take. The main numerical issue con¬ 
cerns the convergence of the Monte Carlo process, which 
is complicated by the exponential proliferation of num¬ 
ber of diagrams with perturbation order a. This leads 
to exponential computational complexity since final re¬ 
sults are recovered only after extrapolation to the infinite 
diagram-order limit. The fermion sign enters the calcula¬ 
tion in an interesting way: different diagrams have differ¬ 
ent signs (arising from the different orderings of fermion 
creation and annihilation operators) and at large diagram 
order the contributions of diagrams with plus and minus 
sign tend to cancel. This cancellation is in fact respon¬ 
sible f or the convergence of the many-body perturbation 
theory.l2lE3I]To manage this issue it is useful to consider a 
Monte Carlo process for approximate series in which the 
maximum diagram order is limited to some finite value 
via a hard or soft cutoff, and then cutoff is increased until 
convergence is reached. 

Diagrammatic Monte Carlo techniques can also take 
advantage of known field-theoretical techniques to run 
the calculation in a self-consistent mode in which cer¬ 
tain infinite series of diagrams are summed and then 
automatically absorbed into the renormalized propaga¬ 
tors and interaction lines using Dyso n-ty pe equations. 
One example is the skeleton expansiorP^ another is a 
“bold” expansion in pertu rbative corrections to an an¬ 
alytic partial resummation.h^^Uas] flexibility allows 

for an additional control over systematic errors coming 
from series extrapolation as well as convergence issues— 
different schemes should produce consistent final results. 

In this work we employ four complementary tech¬ 
niques: 

• A [/-scheme based on a Taylor series expan¬ 

sion of E in powers of U with fixed shifted chemical po¬ 
tential /i = /r — Un/2, see Refs. 1^ and [TMl The total 


electron density, n = n^, and the chemical potential 

are computed posteriori (after results are extrapolated to 
the Umax CO limit). 

• A G^VE-scheme based on skeleton series for E and 
n in which all lines in the diagram are understood as 
fully dressed Green’s functions and screened interactions. 
Self-consistency is implemented by feedback loops when 
G and W are obtained by solving algebraic Dyson equa¬ 
tions, G-i = [G(°)]-i - E and W-^ = [/-^ - R, in 
momentum-frequency representation, see Refs. I135II137I 
for more details. 

• A G^T-scheme based on the skeleton series for E 

and Er when all lines in the graph are understood as 
fully dressed single-particle (Green’s functions) and two- 
particle propagators. This compact formulation is pos¬ 
sible only for a contact interaction potential when the 
sum of ladder-type diagrams for spin-up and spin-down 
particles has the same functional structure as the single¬ 
particle propagator, see Refs. 113211138L and 11391 Again, 
self-consistency is implemented by feedback loops us¬ 
ing Dyson equations, G“^ = — E and T”^ = 

[r(o)]-i - Er, where T^o) is the sum of bare ladder dia¬ 
grams. 

• A [G*'°^]^r(°^-scheme based on diagrams expressed in 
terms of bare single- and pair-propagators with shifted 
chemical potential fl = p — Un/2; this is similar in spirit 
to the [G*^°^]^[/-expansion, but with one extra geomet¬ 
rical series (bare ladder diagrams) being accounted for 
analytically. 

To establish the parameter region where DiagMG 
works, we performed calculations using all four schemes. 
The results were compared to each other and to those 
obtained by DCA. Additional insight was a lso gained by 
doing calculations in the atomic limit We find that 
for bare coupling U/t < 4 and temperature T/t > 0.1, 
all schemes produce consistent results within statistical 
errors. At half filling, n = 1, and U/t < 6, the (G^°))^[/- 
and (G^°^)^r*^°^-schemes still produce results consistent 
with those obtained by DCA and the determinant Monte 
Carlo method. In the dilute region (small filling fac¬ 
tors), the G^T- and (G^°^)^r(°)-schemes can be applied 
for larger values of U. For dilute systems our benchmarks 
include points {U/t = 6,n = 0.6), {U/t = 6,n = 0.3), 
and {U/t = 8,n = 0.3). 


IV. EXTRAPOLATIONS AND 
UNCERTAINTIES 

A. Extrapolations 

All of the numerical methods we have considered rely 
on the extrapolation of results to a thermodynamic or 
asymptotic limit. For DiagMC, which is formulated di¬ 
rectly on an infinite system, the extrapolation is in di¬ 
agrammatic order. All other methods are extrapolated 
from a finite embedded system, finite cluster, or cylin¬ 
der with finite width to the infinite system size limit. 
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FIG. 1. (Color online) Extrapolations of the ground state 
energy at 17/t = 8, n = 1. Main panel: AFQMC and FN 
extrapolated as a function of the inverse cube of the system’s 
linear dimension, L, along with DMRG extrapolated in cube 
of inverse cylinder circumference (also denoted L). DMRG 
data are presented both for rotated (with %/2) and unrotated 
wrapping of cylinders. Inset: DMET data for clusters of size 
and geometry indicated, plotted against the reciprocal of the 
square root of the total number of sites in the cluster N^. 
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In many cases, a considerable contribution to our errors 
comes from this extrapolation procedure, which differs 
from method to method. In some cases determining 
stochastic uncertainties in extrapolation are not possi¬ 
ble, in which case we produce estimates of uncertainties 
by choosing a reference system for a given technique. We 
then assume a normal distribution of uncertainty with re¬ 
spect to the reference. Specifics of the extrapolation pro¬ 
cedure for each system (and of the associated procedure 
for estimating extrapolation uncertainties) are described 
in Sec. m All methods have therefore defined procedures 
to estimate error in TL quantities as accurately as pos¬ 
sible through the use of known reference systems. Ad¬ 
ditional uncertainties due to extrapolation, curve fitting, 
truncation in excitation order are addressed on a per- 
technique basis. These added uncertainties are assumed 
to be normally distributed and defined such that they 
can, in principle, be made arbitrarily small by adding 
additional data. This section illustrates these extrapola¬ 
tions and presents some of the challenges encountered in 
performing them. 

We start our discussion with ground state properties. 
In Fig. [l] DMRG, FN and AFQMC results are presented 
in the main panel. The system size dependence in all 
three techniques is clearly visible and the difference be¬ 
tween the estimated thermodynamic limit value and the 
largest size considered is in some cases outside of the 
error bars of the thermodynamic limit value: in other 
words, extrapolation is essential for obtaining the ther- 
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FIG. 2. (Color online) Energies obtained at U/t = 4, 
T/t = 0.25 and n = 0.8 from G^F (magenta) and 
(turquoise) as a function of the inverse square of the diagram 
order parameter a (upper axis label) along with DCA results 
obtained from finite clusters plotted against the inverse of the 
cluster size (lower axis labels). 


modynamic limit value. For this reason, in the results 
sections we typically present both the thermodynamic 
limit value and the sequence of finite size results which 
led to it, so that the reader can see how large an extrap¬ 
olation is required. 

For DMRG, results from both un-rotated cylinders 
(filled symbol; other smaller-L data are not shown but lie 
on the same scaling curve) and rotated cylinders (crosses) 
are consistent, both scaling as l/-b^, although with dif¬ 
ferent slopes, allowing for a clean extrapolation to the 
thermodynamic limit. On the other hand, the AFQMG 
data indicate a change in scaling for system sizes larger 
than 10 X 10 geometry under twist averaged boundry 
conditions. This could indicate either that unidentified 
complications occur in the large system AFQMC calcu¬ 
lations or that deviations from the 1/L^ size dependence 
might occur in the DMRG data at larger cylinder size 
(i.e. that the DMRG error bar is underestimated). In 
this regard, it is important to note that the ground state 
energy of the largest system examined in DMRG, rotated 
6 X oo, is within uncertainty of the extrapolated AFQMC 
data. 

The FN data also demonstrate a systematic depen¬ 
dence of the energy on system size, allowing a precise 
thermodynamic limit extrapolation. The deviation of 
the FN results from the AFQMC and DMRG results 
is caused by a systematic fixed node error, which by 
comparison to other methods seems to be no more than 
2 X IQ-H. 

Shown in the inset of Fig. are the extrapolations in 
DMET which scale as 1 /y/Nc- Due to the restricted small 
system sizes in DMET and large U, the resulting uncer¬ 
tainty is dominated by the extrapolation. The value, 
also shown in the main panel, is in good agreement with 



















13 


DMRG and FN, and only slightly outside error bars of 
the AFQMC result. 

In Fig. we show data for two DiagMC methods at 
U/t = 4, T/t = 0.25 and n = 0.8, along with DCA 
data. DiagMC is done directly in the thermodynamic 
limit, and the results become successively more precise as 
more and more expansion orders are added to the series. 
The results from the two diagrammatic series we show 
agree within error bars, with the G^F series converging 
more smoothly than The convergence with 

expansion order in the regimes we present is very rapid, 
so that the value at order a = 6 or 7 can be taken as 
representative for the infinite order series, with error bars 
estimated by statistics and by comparison to the results 
at the second largest order; in other words, extrapolation 
to a —oo is not needed. 

DCA for the 2D Hubbard model approaches the ther¬ 
modynamic limit ~ 1/A^c- However, in the parame¬ 
ter regime considered here the many-body physics is 
converged with respect to and deviations from the 
thermodynamic limit are dominated by single-particle 
shell effects. In other regimes, especially at larger 
U, extrapolation in 1/Nc is required, see for example 
Refs. Ea Ell and ETI 

The key result of this section is that in many cases 
extrapolation to the infinite system size limit is needed 
to obtain accurate results, with the value obtained by 
extrapolation significantly different from the value ob¬ 
tained by the largest size studied. For this reason we 
will typically display below both the extrapolated ther¬ 
modynamic limit results and the finite size results that 
produced the extrapolation. 


B. Sources of Uncertainty 

For clarity we repeat the main sources of uncertainties 
and the meaning of the error bars shown in the graphs 
for each technique; further details can be found in the 
sections on each method. 

• AFQMC: at n = I error bars include all sources of 
uncertainty; stochastic errors and extrapolation to TL. 
For n ^ 1, uncertainty from the constrained path ap¬ 
proximation is not estimated by the error bar. 

• FN: error bars account for stochastic Monte Carlo 
errors and for extrapolation to the TL. Uncertainties due 
to the hxed node approximation are not included in the 
error bar. 

• MRPHF: results are not extrapolated to the TL and, 
on each finite system, an estimate of the uncertainty due 
to truncation in the number of Slater determinants is not 
included. 

• UCCSD: error bars do not include an estimate of 
uncertainty for truncation of excitation order to doubles. 

• DMRC: error bars include all sources of uncertainty; 
the extrapolation in the number of basis states and ex¬ 
trapolation to TL. 

• DMET: error bars include all sources of uncertainty; 


uncertainty due to extrapolation in of the number of basis 
states of the impurity solver and extrapolation to TL as 
well as estimates of DMET self consistency convergence. 

• DCA: error bars include all sources of uncertainty; 
stochastic Monte Carlo uncertainties and an additional 
estimate of uncertainty due to extrapolation to the TL. 

• DF: values are presented without error bars, the ef¬ 
fect of neglecting non-ladder and higher order diagrams 
is not quantified. 

• DiagMC: error bars include all sources of uncertainty; 
the stochastic Monte Carlo uncertainty at each expansion 
order and estimate of uncertainty in convergence of ex¬ 
pansion order. 


V. RESULTS AT INTERMEDIATE TO STRONG 
INTERACTION STRENGTH 

We begin our discussion of results with an analy¬ 
sis of an intermediate-to-strongly coupled case, namely 
U/t = 8. Throughout all figures we use common leg¬ 
ends, distinguishing techniques by symbol and color. We 
present both results for the thermodynamic limit and 
the finite system size data from which the thermody¬ 
namic limit results were obtained. This information is 
useful in assessing both the importance of the extrapola¬ 
tion and other aspects of the performance of the method. 
We use open symbols to denote values in the thermody¬ 
namic limit and filled symbols for finite size values from 
which the extrapolations are obtained. 


A. Half-Filled, particle-hole symmetric case 

{U/t = 8, n = 1, t'/t — 0) 

We begin our discussion with an analysis of the energy 
per site for a widely accessible parameter choice, half 
filling, showing in Fig. the temperature dependence of 
the energy and in Fig. fflan expanded view of the T = 0 
energy. 

First we discuss the DCA results in FigH in this 
particle-hole symmetric parameter regime the impurity 
solvers have no sign problem and cubic complexity, mean¬ 
ing that reliable results can be obtained on relatively 
large clusters and that the Monte Carlo errors (here are 
on the order of 10“^) can be systematically reduced with 
additional computation. Thermodynamic limit data are 
obtained from the l/N^. extrapolation. Computational 
scaling towards low temperatures results in an increase 
of uncertainty for fixed computational time, and this is 
reflected in the uncertainty in the extrapolated values. At 
T/t = 0.5 our results agree within error bars with high 
temperature series and lattice Monte Carlo data (see Ref. 

El- 

The results of a DF calculation are shown at T/t = 0.5 
(lower T data are not available). The DF technique ne¬ 
glects vertex functions of higher order than two-particle 
vertices. Further, at the two-particle level we sum only 
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FIG. 3. (Color online) Temperature dependence of the en¬ 
ergy for n = 1 for U/t = 8 obtained by DCA (black circles) 
and DF (red cross) and compared to zero temperature results 
compiled from various techniques. Solid symbols represent 
finite systems, open symbols represent extrapolations to the 
thermodynamic limit (TL). 


a ladder series in the spin and charge channels. Despite 
these approximations, we see that the DF technique pro¬ 
vides an energy which falls on top of that of the extrap¬ 
olated thermodynamic limit DCA result. 

Results from a variety of algorithms are available at 
zero temperature. Fig. presents an expanded view of 
the T = 0 results, with the energy on the vertical axis 
and data for each method offset in the a;-axis. Note that 
in some cases the thermodynamic limit results are further 
offset for clarity. 

We start our discussion of zero-T results with a Monte 
Carlo technique, AFQMC, which is extrapolated to the 
thermodynamic limit. In this case, finite size results are 
averaged over twisted-boundary conditions, which allows 
a smooth and rapid convergence to the thermodynamic 
limit. These results, obtained at half filling from Monte 
Carlo, are unbiased and therefore expected to be exact 
within a quoted uncertainty of ±0.0002t. 

DMRG results on cylinders of infinite length but finite 
width of 3, 4, 5, and 6 for 45 degree rotated systems and 
width of 4 and 6 for non-rotated systems are shown. All 
the finite size data are after phase averaging, showing 
only very weak finite-size effects, so that an extrapola¬ 
tion to the thermodynamic limit is feasible. In this case, 
the estimation of uncertainty (as discussed in Sec. IIIF) 
contains the uncertainty of each extrapolation and the 
difference between the two orientations (rotated and non- 
rotated), both of which are on the order of The 

resulting energy is close to, but slightly outside of, the 
AFQMC results. This extrapolation issue was discussed 
in more detail in Sec. IIV Al 

For DMET we show results obtained for finite clusters 


♦ UCCSD O UCC X AFQMC - TL ► MRPHF 
A DMET - TL ■ DMRG ◄ EN 
A DMET □ DMRG - TL < EN - TL 
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FIG. 4. (Color online) Thermodynamic limit (TL) ground 
state energy for n = 1 for U/t = 8 as obtained by various 
algorithms (open symbols). Also shown are the finite size 
systems (filled symbols with adjacent labels) from which the 
TL ground state energy was obtained. Data from AFQMC 
(red crosses), DMET (blue triangles), UCCSD (maroon dia¬ 
monds), MRPHF (purple triangles), DMRG (orange squares), 
and FN (green triangles). Horizontal thin dotted lines show 
the best estimates for the ground state energy. 


of size 2 X 2, 4 X 2, 8 X 2, and 4x4. The thermodynamic 
limit is obtained by extrapolating the 2 x 2, 4 x 2, and 
4x4 clusters in 1/^/T//,. Errors from the solution of the 
finite impurity are on the order of DMET clus¬ 

ter size convergence is slower at large U, thus U/t = 8 
corresponds to the largest half filling DMET error bar 
discussed here. The total thermodynamic limit uncer¬ 
tainty is estimated to be O.OOlt, and comes entirely from 
the thermodynamic limit extrapolation. The lower end 
of the DMET error bar lies at the average of the DMRG 
and AFQMC estimates. 

For the FN technique, a diffusion Monte Carlo calcula¬ 
tion based on the nodal structure of a trial wave function 
obtained with variational Monte Carlo, we show finite 
size results for a sequence of 45-degree rotated clusters 
with size 98, 162, and 242, which have the property of 
being closed shells at U = 0. The results show only 
weak size dependence, so that the thermodynamic limit 
value shown is close to the finite size results. However, 
the results are systematically above the values obtained 
by AFQMC and DMRG, while they are consistent with 
DMET. This is a consequence of the fixed node approx¬ 
imation, which in this particular case resulted in a fixed 
node error of about 0.0015t. 

The results of UGGSD are shown for systems of size 
6 X 6, 8 X 8, and 10 x 10, and exhibit weak finite size 
effects at this U value. We see that the result is accu¬ 
rate to roughly the 1% level. The deviation is caused by 
correlations that are not captured by singles and dou¬ 
bles. Higher order excitations (triples, quadruples, etc) 
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will eventually recover the remaining energy. To sup¬ 
port this claim we show a single case in Fig. labeled 
as UCCSDT(Q) which includes all triples and a subset 
of quadruples. In this higher order approximation the 
deviation from other techniques is reduced by a factor of 
2. The higher order corrections are more important for 
these coupled cluster calculations than extrapolations in 
cluster geometry size, however, the improvement with 
increased excitation order converges slowly. Also shown 
are approximate results including quadruples from small 
system sizes (4 x 4) which we label UCCSDTQ*. While 
not exact, this approximation scheme produces a results 
which deviates from AFQMC by only 0.15%. 

MRPHF calculations have been performed for several 
finite systems (4x4, 6x6, and 8 x 8) . As summarized 
in the methods description, reaching a constant level of 
accuracy would require a successively larger MR expan¬ 
sion. Results for larger systems are therefore solved less 
precisely; in particular the energy of the 8x8 lattice 
in Fig. 1^ is too high. More sophisticated implementa¬ 
tions and additional optimizations may make it possible 
to reach the accuracy needed to perform extrapolations 
to the thermodynamic limit. 


• DCA A DMET - TL ► MRPHF < FN - TL 
O DCA-TL * DMFT + DF ■ DMRG 
O UCCSD - TL X AFQMC - TL □ DMRG - TL 



FIG. 5. (Color online) Double occupancy data for U/t = 8 
and n = 1. Main panel: temperature dependence of dou¬ 
ble occupancy, obtained from DCA (finite T, black circles) 
and DF (finite T, red plus sign), and the T — 0 techniques 
AFQMC (red crosses), DMET (blue triangles), UCCSD (ma¬ 
roon diamonds), MRPHF (purple triangles), DMRG (orange 
squares), and FN (green triangles). Solid symbols represent 
finite systems, open symbols represent extrapolations to the 
thermodynamic limit (TL). Inset: data at T = 0 reproduced 
with an arbitrary x-axis offset, from MRPHF, UCCSD, FN, 
DMET, DMRG, and AFQMC. 

We now discuss the results for the double occupancy 
in Fig. [^at U/t = 8 and n = I. Open symbols denote 
results in the thermodynamic limit, filled symbols results 
on finite systems. The finite-T DCA results show that 


the double occupancy contribution rises as the tempera¬ 
ture is lowered. The finite-T results are consistent with 
the T = 0 values obtained by AFQMC, DMRG, DMET, 
FN, and MRPHF. At T/t = 0.5, the double occupancy 
obtained from the DF technique is also shown. Unlike 
the total energy, the DF double-occupancy shows devia¬ 
tions from the DCA result, suggesting a cancellation of 
errors in the kinetic and potential energy terms. As for 
all other points we have examined, the DF method pro¬ 
duces results which lie between single-site DMFT values 
(not shown) and the extrapolated DCA results. 

The inset shows the various T = 0 values. Within error 
bars, there is agreement between AFQMC, DMET and 
DMRG results for the double occupancy. DMET obtains 
a value (after thermodynamic limit extrapolation) com¬ 
parable to AFQMC and, overall, shows a weaker system 
size dependence than for the energy. 

UCCSD and FN produce a double occupancy which 
is underestimated as compared to AFQMC and DMET. 
Einite size effects of EN are on the order of O.OOOI. Ei- 
nally, for MRPHE we quote two values for 4 x 4 and 
6x6 systems which show a system size dependence on 
the order of 0.001. This makes a thermodynamic limit 
extrapolation impractical. 


B. Doped strongly correlated regime {U/t = 8, 
n = 0.875, t'/t = 0) 

The half-filled particle-hole symmetric case of Sec. |V A| 
is in many ways ideally suited for numerical algorithms: 
a large charge gap allows methods like the DMRG 
to quickly converge, and particle-hole symmetry makes 
Monte Garlo simulations without a sign problem possi¬ 
ble. We now turn to a case which is particularly difficult 
to simulate, where we expect results to be substantially 
less accurate than for the half filled case. This parame¬ 
ter regime shows metallic behavior, strong particle-hole 
asymmetry, and interesting inhomogeneous phases in the 
ground state. In Pig. we plot the total energy per site 
and in Fig. [^the double occupancy per site at U/t = 8 
and n = 0.875. 

The main panel presents data as a function of temper¬ 
ature. The DCA results remain consistent with T = 0, 
but results are not available at the lowest benchmark 
temperature, T/t = 0.125, due to a large sign problem in 
the Monte Carlo impurity solver. The finite size effects, 
at the system sizes accessible in DCA, are smaller than 
in the n = 1 case. 

The inset to Fig. [^presents data at T = 0 with an arbi¬ 
trary cc-axis offset added for clarity. The AFQMC simula¬ 
tion, using a (non-variational) constrained path approx¬ 
imation in the absence of particle-hole symmetry, yields 
a result for the total energy that is lower than the one 
obtained from DMET, PN, and DMRG. The total energy 
difference is ~1% when compared to finite sized DMRG, 
and ^ 1.4%(2.1%) when compared to the DMET 8x2 
(thermodynamic limit) cluster respectively, and ^ 2.1% 
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FIG. 6. (Color online) Data for n — 0.875 for U/t = 8. Main 
panel: temperatnre dependence oi Ejt compiled from vari¬ 
ous techniques. Solid symbols represent finite systems, open 
symbols represent extrapolations to the thermodynamic limit 
(TL). Finite T resnlts are shown for DCA (black circles), and 
zero-T data from AFQMC (red crosses), DMET (blue trian¬ 
gles), UCCSD (maroon diamonds), DMRG (orange squares), 
and FN (green triangles). Top left inset: zoom in to the 
zero-T data from DMET, AEQMC, FN, and DMRG includ¬ 
ing finite system size data (as labeled). 

in comparison to FN. 

DMRG shows the results for cylinders of infinite length 
and finite widths of 4 and 6 lattice sites after using phase 
averaging. The energy is higher for the wider cylinder, 
and for the width 6 cylinder the energy is above the en¬ 
ergy from AFQMC. Given that the extrapolation is per¬ 
formed with only two widths, we consider the extrapo¬ 
lated DMRG value to be not reliable in this case and 
omit it entirely. 

DMET shows a large system size dependence and a de¬ 
pendence of the thermodynamic limit value on the cluster 
sequence chosen for the extrapolation. We show an ex¬ 
trapolation based on 2 x 2, 4 x 2, and 8x2 clusters. The 
use of the 8x2 cluster allows inhomogeneous order to de¬ 
velop, giving an extrapolated value oi Ejt = —0.749(7). 
The extrapolation using the 4x4 rather than 8x2 cluster, 
which does not allow for inhomogeneous order, yields a 
value oi Ejt = —0.737(5). Since the energy changes non- 
monotonically: the 8 x 2 energy lies above the 2x2 energy, 
but below the 4x2 energy, the uncertainty in the thermo¬ 
dynamic limit extrapolation is very large, and does not 
provide any more information than the results obtained 
from the largest clusters. 

The FN method shows a clear finite system size depen¬ 
dence. The infinite system value is estimated from the 
16 X 16 and 20 x 20 values, and finite size errors are on 
the order of O.OOK, much larger than the stochastic er¬ 
rors of O.OOOOlt. Here, the FN results are consistent with 
DMET extrapolation which omits the 4x4 cluster and 


these are considerably higher than AEQMC. This is sug¬ 
gestive of a fixed node error of « 0.015t, indicating that a 
uniform variational wave function may not be enough to 
fully account for the nature of the ground state. Indeed, 
the VMC error is of the order of 0.022t, much larger than 
the one obtained at half filling which was « 0.004t. In 
both cases, the EN projection improves the VMC results 
by the same order of magnitude. 

UCCSD and MRPHE results are much higher in en¬ 
ergy {E/t = —0.7094(5) for MRPHE (16 x 4 system) 
and E/t = —0.7122 for UCCSD, barely visible on the 
main panel), an indication that correlated metallic states 
are difficult to capture with these methods. Although 
not shown, data is available for UCCSD(T) (perturba¬ 
tive inclusion of triples) in the supplemental material,^^ 
which improves upon the value from UCCSD, and gives 
E/t = -0.7272 (-0.7281) for a 16 x 4 (16 x 8) cluster. 
Full inclusion of triples (UCCSDT) lowers the 16 x 4 es¬ 
timate to E/t = —0.7427. The MRPHF results indicate 
the need for a much larger MR expansion than that af¬ 
forded in this work. 

In this parameter regime, ordered ‘stripe’ phases might 
exist. However, the precise form of these stripes is 
strongly influenced by choice of finite size systems (e.g. 
width and orientation of the cylinder in DMRG and 
shape of the cluster in DMET) that are used for the ther¬ 
modynamic extrapolation and the approximations used 
to solve that finite system. The finite temperature al¬ 
gorithms have not reached the onset of inhomogeneous 
states at the lowest temperature accessible. The precise 
nature of the inhomogeneities in the ground state in this 
parameter regime is still open. 

We finally briefly mention the results for double occu¬ 
pancy in Fig. for U/t = 8 and n = 0.875. As was the 
case for the energies, the finite-T results smoothly con¬ 
nect to the zero-T values. MRPHF and UCCSD overesti¬ 
mate the double occupancy by close to 15%. The remain¬ 
ing ground state methods (DMRG, AFQMC, DMET, 
and FN) present consistent values in the range between 
0.04 and 0.043. Both FN and AFQMC values contain 
additional (fixed node and constrained path) errors that 
are not estimated by the error bar. 


C. Half-filled, non particle-hole symmetric case 

{U/t = 8, n = 'i-,t!/t = —0.2) 

We now turn our attention to a case of half filling with¬ 
out particle-hole symmetry, by adding a second nearest 
neighbor hopping t'. An overview of the energies from 
several algorithms for U/t = 8, n = 1 and t' = —0.2 is 
shown in Fig. 

The main panel shows the temperature-dependence of 
the data. The DCA results available at finite T show 
almost no sign problem for T/t — 0.5 and T/t = 0.25, but 
are hampered by a severe sign problem at T/t = 0.125. 
The results are consistent within error bars with the zero- 
temperature results. 
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FIG. 7. (Color online) Data for n = 0.875 for U/t = 8. 
Main panel: temperature dependence of double occupancy, 
D, compiled from various techniques. Solid symbols repre¬ 
sent finite systems, open symbols represent extrapolations to 
the thermodynamic limit (XL). Finite T resnlts are shown 
for DCA (black circles), and zero-T data from DMET (blue 
triangles), MRPHF (purple triangles), UCCSD (maroon dia¬ 
monds), DMRG (orange squares), and FN (green triangles). 
Inset: zoom in to the zero-T data from FN, AFQMC, DMET, 
and DMRG. 


As at I7/t = 8, n = 0.875, t' = 0, the AFQMC is ap¬ 
proximate because of a constrained path approximation 
due to the lack of particle-hole symmetry. Despite this, 
the results are in agreement with both the DMET and 
DMRG results. 

The DMET results are obtained on clusters of size 
2 X 2, 2 X 4, and 4x4. Errors of the individual finite 
size systems are substantially smaller than the system 
size dependence. The DMET thermodynamic limit is 
consistent with the thermodynamic estimates obtained 
from DMRG (from cylinders of width 4 and 6) and from 
AFQMG. This is even more evident in the bottom right 
inset, which displays the thermodynamic limit estimates 
on a smaller scale. 

FN results are higher in energy than AFQMC and 
DMET (well within two joint standard deviations) and 
are higher than DMRG by O.OOlSt. As seen in previous 
plots, the finite system size dependence of the fixed node 
results is small on this scale. 

UCCSD results show only small finite size effects and 
an overall energy «1% higher than other techniques. The 
MRPHF results obtained on finite systems show an en¬ 
ergy that rises rapidly as the system size is increased. As 
in the case of = 0, a systematic extrapolation to the 
thermodynamic limit is not possible, and we only present 
results on finite systems. 


♦ UCCSD X AFQMC - XL ► MRPHF 

* DMEX •< FN ■ DMRG • DCA 

A DMEX-XL < FN-XL □ DMRG-XL O DCA - XL 



FIG. 8. (Color online) Data for n = 1 for U/t = 8 with t' jt = 
—0.2. Main panel: temperature dependence of E/t compiled 
from various techniques. Solid symbols represent finite sys¬ 
tems, open symbols represent extrapolations to the thermody¬ 
namic limit (TL). Finite T results are shown for DCA (black 
circles), and zero-T data from AFQMC (red crosses), DMET 
(blue triangles), UCCSD (maroon diamonds), MRPHF (pur¬ 
ple triangles), DMRG (orange squares), and FN (green trian¬ 
gles). Top left inset: zoom in to the zero-T data from MR¬ 
PHF, UCCSD, DMET, FN, DMRG, and AFQMC, including 
finite system size data (as labeled) for MRPHF, FN, DMET, 
and DMRG. Bottom right panel: enlarged region of the top 
left inset showing DMET, DMRG, FN, and AFQMC data 
at T = 0, including error bars, extrapolated to the infinite 
system size. 


VI. RESULTS IN THE INTERMEDIATE 
COUPLING REGIME 

In this section we repeat the previous discussion for an 
interaction strength of half the size, U/t = i. As before 
we start our discussion at half filling. We then discuss a 
correlated metallic case with 20% doping. 

A. Half-Filled, particle-hole symmetric case 

(U/t = 4, n = 1, t!/t — 0) 

In Fig. we report the energy as a function of temper¬ 
ature. At finite T and U/t = 4, both DGA and the di¬ 
agrammatic Monte Carlo method for the series 

provide results in the thermodynamic limit. DCA results 
in the thermodynamic limit are extrapolated from finite 
clusters, DiagMC results are extrapolated in the expan¬ 
sion order. The results are consistent within the error 
bars of the respective methods. The large error bars of 
the extrapolation in DiagMC-mainly come from 
a conservative estimate of the diagram-order extrapola¬ 
tion error. DCA shows surprisingly large finite size ef¬ 
fects which persist above Nc = 72, unlike at U/t = 8. 
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♦ UCCSD >• MRPHF 
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FIG. 9. (Color online) Data for n = 1 for U/t = 4. Main 
panel: temperature dependence oi Ejt compiled from vari¬ 
ous techniques. Solid symbols represent finite systems, open 
symbols represent extrapolations to the thermodynamic limit 
(XL). Finite T results are shown for DCA (black circles) and 
DiagMC (turquoise stars), and zero-T data from AFQMC 
(red crosses), DMET (blue triangles), UCCSD (maroon dia¬ 
monds), MRPHF (purple triangles), DMRG (orange squares), 
and FN (green triangles). Top left inset: zoom in to the zero- 
T data from MRPHF, DMET, FN, DMRG, and AFQMC, 
including finite system size data (as labeled) for MRPHF, 
FN, DMRG, and DMET. 


FIG. 10. Data for n = 1 for U/t = 4. Main panel: tem¬ 
perature dependence of double occupancy, D, compiled from 
various techniques. Solid symbols represent finite systems, 
open symbols represent extrapolations to the thermodynamic 
limit (TL). Finite T results are shown for DCA and DiagMC 
(turquoise stars), and zero-T data from DMET (blue trian¬ 
gles), UCCSD (maroon diamonds), MRPHF (purple trian¬ 
gles), DMRG (orange squares), and FN (green triangles). In¬ 
set: zoom in to the zero-T data from DMRG, FN, UCCSD 
(only finite systems data), DMET, and AFQMC. 


While each individual result has uncertainties in the 
energy on the order of the spread in values results 

in large uncertainty when extrapolated to the thermody¬ 
namic limit. 

We now move to the zero-temperature methods, which 
are shown in the inset of Fig. AFQMC provides numer¬ 
ically exact ground state energies for this system. The 
value quoted is E/t = —0.8603(2), which is in agreement 
with the DMET value of E/t = —0.8604(3). and the 
DMRG value of E/t = —0.8605(5). DMET values are 
obtained from an extrapolation of 2 x 2, 2 x 4, and 4x4 
clusters. DMRG values are obtained from an extrapola¬ 
tion of widths 3, 4 and 5 for 45-degree rotated cylinders 
and of width 4 and 6 for non-rotated cylinders. 

The results obtained by AFQMC, DMRG, and DMET 
are in excellent agreement with recent calculations ob¬ 
tained from linearized auxiliary fie lds Monte Carlo 
(LAQMC) available in the literatur^SD^ which gives 
E/t = —0.85996(5). FN results are higher in energy 
{E/t = —0.8575(3)), and unlike in previous cases for 
stronger interaction, a clear dependence on the finite sys¬ 
tem studied is visible. The FN projection technique leads 
to an energy gain of « 0.002t with respect to the VMC re¬ 
sult oi E/t = —0.8558(5). This number can be compared 
with a previous estimation of the thermodynamic limit 
in VMC obtained with a slightly less accurate variational 


state, see Ref. Ei 


The UCCSD thermodynamic limit overestimates the 
energy by ^ 0.7%. MRPHF values, shown as purple right 
triangles in the main panel and inset, show large finite 
size effects and are higher than the values obtained with 
other methods. We see that as system size is increased, 
the energy increases rapidly. 


In Fig. we report the double occupancy vs T. At 
finite T, the DCA results show a clear rise in £> as T de¬ 
creases from 0.5t —>■ 0.25L However, this trend reverses 
as T decreases further. A similar behavior is obtained 
also by using DiagMC, demonstrating that this is a gen¬ 
uine effect present in the Hubbard model. 


These trends are consistent with the T — 0 data, which 
lie below all of the DCA data points at finite T. For 
clarity of presentation we again display this data in the 
inset and add an arbitrary x-axis offset. We see that fi¬ 
nite size effects in DMET are very small, and that the 
extrapolation of DMET agrees perfectly with AFQMC. 
Finite sized FN results produce values comparable to to 
DMRG, DMET and AFQMG. However, extrapolation in 
FN results in an underestimate of the double occupancy, 
although within uncertainty. In the DMRG simulations, 
phase averaging has greatly reduced finite size effects, 
and the DMRG error bars are determined by the trunca¬ 
tion errors. Within those error bars, DMRG results are 
consistent with AFQMC and DMET. 
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B. Doped case {U/t = 4:,n = 0.8, t'jt = 0) 


Away from half filling (with t' = 0), we can perform 
further comparisons at finite T between DCA and Di- 
agMC at U/t = A. We begin the discussion with the inset 
of Fig. [m which shows the convergence of the imaginary 
part of the local Matsubara self-energy of the G^F and 
G^W DiagMC series as a function of evaluation order. 
The values are compared to DCA results. We see that 
the first six orders of the series are precise enough to 
get good agreement of the Matsubara self-energy in the 
thermodynamic limit, and convergence is rapid. While 
deviations are visible in the energy, these are attributed 
to differences of the chemical potential, i. e. the real part 
of the self-energy. 


The top left inset of Fig. 11 shows the convergence of 
the energy in DCA and two different series of DiagMC, 
[( 7 (o)] 2 fy G^F, for increasing order of the diagram¬ 
matic resummations, a. The values obtained from the 
three techniques are within error bars. 

The T > 0 values smoothly connect to T = 0 (although 
a precise comparison is not possible because we lack a 
quantitative functional form to extrapolate the T > 0 
values to T = 0) which we display separately in Fig. 
where data from DMET, AFQMC (constrained path]7 
UCCSD, MRPHF, FN, and DMRG are shown. In this 
case, MRPHF and UCCSD are systematically higher in 
energy from the other techniques. In the case of UCCSD 
we see larger finite size effects than at U/t = 8. Inclu¬ 
sion of higher orders of excitation (perturbative triples 
(T), triples T and perturbative quadruples (Q)) suggests 
that the dominant error is associated with the truncation 
of the excitation order and not finite size effects. The 
FN result is in agreement with the value from DMRG. 
At slightly lower energy AFQMC (constrained path) and 
DMET are in close agreement. Overall, the spread in 
energies is similar to that shown away from half filling 
at U/t = 8 (see Fig. but smaller in magnitude, per¬ 
haps due to better convergence for weak coupling in some 
techniques. 


VII. RESULTS IN THE WEAK COUPLING 
REGIME 


In the weak coupling limit we restrict the presentation 
of data to the half-filled case since the correlated metallic 
phase is not qualitatively distinct from U/t = 4. Data 
sets for doped, weakly correlated systems are available in 
the supplemental material.!® 


A. Half-Filled, particle-hole symmetric case 

(U/t = 2, n = 1, t' jt — 0) 


In Figs. 13 and 14 we present results for U/t = 2 and 
n = 1, the half-filled weak coupling regime. This regime 
is particularly easy for methods based on a weak coupling 


* DiagMC gT finite order * DiagMC G'F inf order 
DiagMC [G"’’] V'" A DMET - TL X AFQMC - TL 

• DCA A fN - TL □ DMRG - TL 

O DCA-TL ► MRPHF ♦ UCCSD 



FIG. 11. (Color online) Data for n = 0.8 for U/t = 4. Main 
panel: temperature dependence ot E/t compiled from vari¬ 
ous techniques. Solid symbols represent finite systems, open 
symbols represent extrapolations to the thermodynamic limit 
(TL). Finite T results are shown for DCA (black circles) and 
DiagMC (pink and turquoise asterisks), and zero-T data from 
AFQMC (red crosses), DMRG (orange squares), FN (green 
triangles), and DMET (blue triangles). Top left inset: zoom 
in to the T/t = 0.25 data from DCA and two types of Di¬ 
agMC for different orders a = 3,4,.... Bottom right inset: 
plot of the imaginary part of the local self energy ImE(ia;„) 
from DiagMC for different expansion orders a and from DCA 
(black circles covered by magenta stars). 
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FIG. 12. (Color online) Thermodynamic limit (TL) ground 
state energy for n = 0.8 for U/t = 4 as obtained by vari¬ 
ous algorithms (open symbols). Also shown are the finite size 
systems (filled symbols with adjacent labels) from which the 
thermodynamic limit ground state energy was obtained. Data 
from AFQMC (red crosses), DMET (blue triangles), UCCSD 
(maroon diamonds), UCC on a 10 x 4 lattice (shaded dia¬ 
monds), MRPHF (purple triangles), DMRG (orange squares), 
and FN (green triangles). Horizontal thin dotted lines show 
the best estimates for the ground state energy. 
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DiagMC [G'"f U O UCCSD - TL O DCA - TL 

XAFQMC-TL bDIVIRG < FN - TL ADMET-TL 
+ DF ■<FN "DCA ■‘DMET ♦UCCSD ►MRPHF dDMRG 



much larger finite size effects, approaching the thermody¬ 
namic limit energy from below (only the largest system 
size is visible on the scale of the plot). DMRG results 
with phase averaging are precise even at U/t = 2, though 
much larger uncertainties than at U/t = 8 are present. 

The results of MRPHF, outside of the scale shown by 
the inset, show a gradual decrease of the energy with in¬ 
creasing system size: 4x4 yields E/t = —1.1260, 6 x6 
E/t = -1.1515, and 8 X 8 E/t = -1.1629. UCCSD re¬ 
sults, here shown for cluster sizes of 12 x 12, 10 x 10, 
8x8, and 6x6, show large finite size effects. With the 
aid of an extrapolation the value in the thermodynamic 
limit is estimated higher than other techniques. The de¬ 
viation of the thermodynamic limit value from AFQMC 
and DMET is on the order of 2 x 10“^t, suggesting that 
excitations beyond the singles and doubles level are im¬ 
portant even in this relatively weak coupling regime. 


FIG. 13. (Color online) Data for n = 1 for U/t = 2. 
Main panel: temperature dependence of E/t compiled from 
various techniques. Solid symbols represent finite systems, 
open symbols represent extrapolations to the thermodynamic 
limit (TL). Finite T results are shown for DCA (black circles) 
and DiagMC (blue asterisks), and zeio-T data from AFQMC 
(red crosses), DMET (blue triangles), UCCSD (maroon dia¬ 
monds), MRPHF (purple triangles), DMRG (orange squares), 
and FN (green triangles). Top left inset: zoom in to the zero- 
T data from UCCSD, DMET, FN, DMRG, and AFQMC, and 
DMET. Bottom right panel: enlarged region of the top left 
inset showing thermodynamic limit data for DMRG, DMET, 
UCCSD, FN and AFQMC data at T = 0, including error 
bars, extrapolated to the infinite system size. 


expansion around a non-interacting system, and many of 
the algorithms show uncertainties that are much smaller 
than in the intermediate or strong interaction limit. 

At non-zero T, the data from two types of DiagMC and 
from DCA in the thermodynamic limit agree within un¬ 
certainty. The values smoothly connect to the T = 0 val¬ 
ues, except for MRPHF energies, which are higher than 
the ground state energies obtained by the other meth¬ 
ods and higher than the energies obtained for the lowest 
temperature point obtained from both finite-T methods. 

At T/t = 0.5 we show a result from DF. As was the 
case in the strong coupling regime, the DF procedure pro¬ 
duces an energy estimate consistent with DCA results. In 
this case, with only weak finite size dependence, the un¬ 
derlying DMFT approximation differs from DCA by only 
0.4%. The DF value improves on the DMFT and differs 
from the extrapolated DCA results by only 0.07%. 

In the lower right inset of Fig.[^ we present T = 0 ex¬ 
trapolations. As in the case of larger interactions, DMET 
and AFQMC (which is numerically exact in this situa¬ 
tion) agree precisely while FN is slightly higher in energy 
but compatible within two error bars. In the upper left 
inset we explore the finite size effects of the methods. 
In the case of DMET, these finite size effects are small, 
and can be extrapolated with small error bars. FN shows 
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FIG. 14. Data for n=l for U/t = 2. Main panel: tempera¬ 
ture dependence of double occupancy, D, compiled from var¬ 
ious techniques. Solid symbols represent finite systems, open 
symbols represent extrapolations to the thermodynamic limit 
(TL). Finite T results are shown for DCA (black circles) and 
DiagMC (turquoise asterisks), and zero-T data from DMET 
(blue triangles), UCCSD (maroon diamonds), MRPHF (pur¬ 
ple triangles), DMRG (orange squares), and FN (green trian¬ 
gles). Inset: zoom in to the zero-T data from FN, UCCSD, 
MRPHF, DRMG, DMET, and AFQMC. 


Finally, Fig. [T^ shows the double occupancy for these 
parameters. We see an increase in finite size effects in 
DCA as we progress to lower temperature. Reasonable 
agreement with DiagMC is achieved in the double occu¬ 
pancy. At temperatures lower than our lowest tempera¬ 
ture, the double occupancy will need to dip, as was the 
case at C/// = 4, in order for the finite T data to be con¬ 
sistent with T = 0. Similar to the strong coupling case, 
DF provides only a slight shift to the double occupancy, 
a minimal improvement over DMFT alone. 

The ground-state double occupancies are very precisely 
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FIG. 15. (Color online) Imaginary part of the local self 
energy, ImE;oc(*a;„) at U/t = 2, Tjt = 0.5 and n = 0.3 from 
DCA and DiagMC. In the case of DiagMC-G^F we label a, 
the series order from Eqn Q. 


determined by AFQMC and DMET which are in agree¬ 
ment. The FN value is somewhat overestimated. Results 
from DMRG fall below AFQMC and DMET and the er¬ 
ror bar underestimates the uncertainty. The larger er¬ 
ror appears consistent with the difficulty in treating the 
small U limit in the DMRG calculations. The results 
from MRPHF show an improvement in D as the system 
size is increased, consistent with the behavior for the en¬ 
ergy. In the case of UCCSD, since it is an expansion in 
the coupling strength, at weak coupling the procedure is 
more reliable, and while there are substantial finite size 
effects, the extrapolation produces a result within error 
bars of AFQMC, and in general agreement with DMET. 


VIII. FREQUENCY AND MOMENTUM 
DEPENDENCE 


Next we discuss single-particle finite temperature prop¬ 
erties. All finite-temperature algorithms discussed in this 
work are based on approximations of the single-particle 
self-energy. We show three characteristic plots for this 
quantity: Figure [Ts] shows the imaginary part of the lo¬ 
cal self-energy as a function of Matsubara frequency, Fig- 
ure[T6| shows the dependence of the real part of the lowest 
Matsubara frequency on k-space, and Figure[^shows the 
frequency-dependence of the imaginary part of the self¬ 
energy for a specific momentum. Any discrepancy in the 
energy or double occupancy is the result of discrepancies 
in the single-particle self-energy. 

The data shown in Fig. [^is obtained for weak inter¬ 
action strength U/t = 2 and for a density n = 0.3. In 
this metallic regime, self-energies are small. Black circles 
denote the imaginary part of the local self energy from 
an Nf. = 20 DCA calculation, which for these parameters 
shows essentially no finite size effects. The data agrees 


FIG. 16. (Golor online) Gomparison of the lowest Matsubara 
frequency, ioJo, real part of the self energy, ReE(fc,ia;o), ob¬ 
tained from DF (red) compared to 72-site DGA calculations 
(black) plotted as a function of momentum, k, throughout 
the Brillouin zone for n= 1.0, U/t = 8, Tjt = 0.5. The Dual 
Fermion and DCA self-energies are plotted as step functions. 
Also included are interpolated r esults obtaine d by diagram¬ 
matic determinantal Monte (dashed black) 

with a gray shading to indicate the level of uncertainty. 


perfectly with DiagMC-G^FF data shown as red dashed 
lines, and convergence of the DiagMC-G^F method to the 
result of the other two methods (stars, magenta dotted 
line) is observed as a function of expansion order a. This 
agreement implies that the local physics is captured well 
by all three algorithms. 


In Fig.j^we examine momentum dependent data. We 
show a path of (kx,ky) through the Brillouin zone and 
plot the real part of the self energy at the lowest Mat¬ 
subara frequency atl7/t = 8, /3 = 2 and n = 1. Di¬ 
agMC data is not available in this regime, but for com¬ 
parison we plot large DCA cluster results (W = 72) 
and results from continuous-time lattice Monte Carlo 
simulations (DDMC, see Refs. 112011140111421 and I143II . 
The DCA approximation (blue lines) produces a step- 
discretized self-energy which is in approximate agreement 
with the momentum dependence from other techniques. 
Discrepancies between the approximate DF method and 
the (essentially converged) DDMC method are visible 
but within the uncertainty of the DDMC comparison 
data. Any discrepancies are expected to rapidly dis¬ 
appear at higher Matsubara frequencies, which can be 
seen in a comparison of the imaginary part of the self¬ 
energy (Fig. 17 1 of DF and DCA at fixed momentum k = 
(tt, 0). For comparison/verification purposes, we include 
results from dynamical vertex approximation (DTA, see 
Refs. 1124111441 and [145]) which, in a spirit similar to DF 
solves the model in an expansion of two-particle vertex 
functions. We find that the results from DF and DTA 
are consistent. 
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FIG. 17. (Color online) Comparison of the frequency depen¬ 
dence of the imaginary part of the self energy, ImE(fc,icu„), 
at fixed k = (tt, 0) obtained from DF (red) compared to 72- 
site DCA calculations (black) plotted for n = 1.0, U/t = 8, 
T/t = 0.5. Also shown ar e result s from dynamical vertex 
approximation, DTA, (blue). ^^ * ^^'^ ^ 


IX. TABLES FOR GROUND STATE 
PROPERTIES 

We conclude the discussion of our results with a list 
of the thermodynamic limit estimates for the half-hlled 
t' = 0 case. Table |n] shows a list of energies in the 
thermodynamic limit, obtained from AFQMC, DMRG, 
DMET, and FN. The MRPHF and UCCSD values pre¬ 
sented show the value for the largest finite size system 
studied. The uncertainties presented are the best uncer¬ 
tainties available within each algorithm and do not con¬ 
tain an assessment of systematic errors (e.g. fixed node 
or truncation in expansion order). We see that for much 
of phase space, errors are a few times 10“‘^t and values 
between the techniques are remarkably consistent. 

Table |III| shows the double occupancy for the same 
values. Relative errors for the double occupancy are of 
the same magnitude as for the total energy, and values in 
the thermodynamic limit are consistent within error bars. 
Tables |IV] and 0 present ground state energies for the 
densities n = 0.8 and 0.875 for values oiUft considered in 
this work. The full table of values for the data presented 
in this paper is available online as supplemental material 
to this paper.l^ 

Not all quantities are as consistent as the energies. 
This is especially true for order parameters and correla¬ 
tion functions, where discrepancies outside of error bars 
are present. Presumably, many competing phases exist in 
a narrow energy window near the ground state, and the 
most favorable state found in each method will depend on 
details of the finite size system and the approximation. 
Table |VI| shows the comparison between the magnetiza¬ 
tion that DMET observes and the magnetization found 
in AFQMC for the full range of U/t at half filling. At 


weak interaction strength, DMET finds a larger polar¬ 
ization than AFQMC even though the energies agree to 
all significant digits, as a result of DMET finite size scal¬ 
ing from small clusters. At large interaction strength, 
AFQMC gives a polarization with large statistical fluc¬ 
tuation despite very accurate energies. Similar behavior 
(not shown here) is apparent for other variables, e.g. the 
d-wave order parameter or the stripe geometry observed 
at U/t = 8 and n = 0.875. 


X. CONCLUSIONS 


In this paper we have presented a detailed examination 
of results for static and dynamic properties of the two di¬ 
mensional Hubbard model at correlation strengths rang¬ 
ing from weak to intermediate to strong coupling, and 
at various carrier concentrations, obtained using state of 
the art numerical methods. We believe the results are 
useful for two reasons. First, the two dimensional Hub¬ 
bard model is one of the paradigm models of quantum 
condensed matter theory, and it is therefore important 
to determine, as reliably as possible, the state of our 
knowledge about it. Second, solving the grand-challenge 
problem of determining the physics of interacting many- 
electron systems will require numerics, and as no one 
technique is likely to provide solutions in all regimes or 
for all quantities of physical interest, it is important to 
develop tools for assessing the strengths and weaknesses 
of different approaches. 

We argue that the only quantities that can meaning¬ 
fully be compared between different approaches are esti¬ 
mates, with error bars, for the thermodynamic limit val¬ 
ues of observables including local operators such as the 
energy, double occupancy, density (or chemical poten¬ 
tial) and magnetization, as well as correlation functions 
such as the electron self-energy. We restricted attention 
to methods and regimes for which large enough systems 
can be studied that reasonable extrapolations to the ther¬ 
modynamic limit can be performed. Care is required in 
performing the extrapolations and we have found it use¬ 
ful to present both the extrapolated results and (in most 
cases) the finite size data that led to the extrapolation. 

Comparison of results obtained from different methods 
shows that the ground state properties of a substantial 
part of the Hubbard model phase space are now under 
numerical control (see e.g. Figs and 12). Moreover, 
where there is agreement on the ground state properties, 
the non-zero temperature methods appear to connect 
smoothly to the ground state as the temperature is de¬ 
creased, although a quantitative extrapolation to T = 0 
is not yet available. The most substantial uncertainties 
exist at intermediate correlations (e.g. U/t « 4 —>• 8) 
and at dopings near to but not equal to the half hlling 
value n = 1. In this intermediate coupling/near half 
filling regime several physically different states seem to 
have very similar energies, and small effects can favor the 
choice of one state over the other, leading to substantial 
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u 1 

2 

4 

6 

8 

12 

AFQMC 

-1.1763 0.0002 

-0.8603 

0.0002 

-0.6568 

0.0003 

-0.5247 

0.0002 

-0.3693 

0.0002 

DMET 

-1.1764 0.0003 

-0.8604 

0.0003 

-0.6562 

0.0005 

-0.5234 

0.0010 

-0.3685 

0.0010 

DMRC 

-1.176 0.001 

-0.8605 

0.0005 

-0.6565 

0.0001 

-0.5241 

0.0001 

-0.3689 

0.0001 

FN 

-1.175 0.001 

-0.8575 

0.0003 

-0.6551 

0.0001 

-0.52315 

0.00005 

-0.36835 

0.00005 

MRPHF 

-1.1628 [8 X 8] 

-0.8554 

[8 X 8] 

-0.6512 

[8 X 8] 

-0.5169 

[8 X 8] 

-0.3626 

[8 X 8] 

UCCSD 

-1.1735 0.0004 

-0.8546 

[14 X 14] 

-0.6510 [10 X 10] 

-0.5191 

[10 X 10] 

-0.3647 

[10 X 10] 

UCCSDTQ* 

-1.1749 - 

-0.8610 

- 

-0.6582 

- 

-0.5255 

- 

-0.3696 

- 


TABLE II. Zero-temperature energy and uncertainty for n = 1, T = 0, for a range of interaction strengths U, obtained from 
AFQMC, DMET, DMRG, FN, MRPHF and UCCSD. Where extrapolations to the TL are not available, finite size geometries 
are listed in lieu of uncertainties. UCCSDTQ* data estimates higher order corrections by including triples from a [6 x 6] and 
quadruples from a [4 x 4]. UCCSD data for [//t > 4 provides nearly converged energy estimates with respect to system size. 


U 

2 

4 

6 

8 

12 

AFQMC 

0.1923 0.0003 

0.1262 

0.0002 

0.0810 

0.0001 

0.0540 

0.0001 

0.0278 

0.0001 

DMET 

0.1913 0.0004 

0.1261 

0.0001 

0.08095 

0.00004 

0.05398 

0.00007 

0.02780 

0.00003 

DMRC 

0.188 0.001 

0.126 

0.001 

0.0809 

0.0003 

0.0539 

0.0001 

0.0278 

0.0001 

FN 

0.198 0.001 

0.125 

0.001 

0.0803 

0.0002 

0.0535 

0.0001 

0.0278 

0.0002 

MRPHF 

0.1824 [8 X 8] 

0.1262 

[8 X 8] 

0.0818 

[8 X 8] 

0.0544 

[8 X 8] 

0.0275 

[8 X 8] 

UCCSD 

0.194 0.002 

0.1268 [12 X 12] 

0.0807 

[10 X 10] 

0.0537 

[10 X 10] 

0.0276 

[10 X 10] 


TABLE III. Zero-temperature double occupancy and uncertainty for n = 1, for a range of interaction strengths U, obtained 
from AFQMC, DMET, DMRC, FN, MRPHF, and UCCSD. Where extrapolations to the TL are not available, finite size 
geometries are listed in lieu of uncertainties. 


uncertainties in physical quantities. Also, the issues asso¬ 
ciated with fermion sign problems seem to be most severe. 
Interestingly, it is this regime that is of most physical in¬ 
terest in connection with high-Tc superconductivity in 
the copper-oxide materials. 

Where two or more methods produce results that agree 
within reasonable error bars, we take the result to be es¬ 
tablished, and appropriate for use as a benchmark. Ta¬ 
bles of our benchmark results are contained in the supple¬ 
mental materiaP^ and made available online. We expect 
that these results will be useful in validating new meth¬ 
ods, or new implementations of existing methods. 

Turning now to prospects and open questions, we first 
observe that all of the methods we have considered to 
date have difficulty in the physically interesting inter¬ 
mediate coupling, near half filling regime. Development 
of new methods, or improvement of existing methods to 
deal with this regime is urgently needed. Further, we 
remark that our understanding of dynamical correlation 
functions, even ones as simple as the electron Green func¬ 
tion, is much less advanced than our understanding of 
ground state properties and simple expectation values. 
Finally, we observe that the process of producing this pa¬ 
per, in particular the confrontation of each method with 
the body of related work produced by other methods, 


led in several cases to substantial improvements in algo¬ 
rithm and error analysis. We suggest that as quantitative 
numerics in the many-electron field continues to evolve, 
intercomparison of different methods, leading to bench¬ 
marking on important model problems, will significantly 
advance the field. 
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u 

2 

4 

6 

8 

AFQMC 

-1.306 

0.002 

-1.110 

0.003 

- 

- 

- 

DMET 

-1.3062 

0.0004 

-1.108 

0.002 

-0.977 

0.004 

-0.88 

0.03 

DMRG 

- 

- 

-1.104 

0.0014 

- 

- 

- 

- 

FN 

-1.3044 

0.0007 

-1.1032 

0.0007 

-0.967 

0.001 

-0.877 

0.001 

MRPHF 

-1.2931 

[10 X 4] 

-1.0910 

[10 X 4] 

-0.9454 

[10 X 4] 

-0.8415 

[10 X 4] 

UCCSD 

-1.3065 [10 X 16] 

-1.0868 

[10 X 16] 

-0.9300 [10 X 16] 

-0.8233 [10 X 16] 

UCCSDT* 

-1.3078 

[10 X 4] 

-1.0981 

[10 X 4] 

-0.9607 

[10 X 4] 

-0.8641 

[10 X 4] 


TABLE IV. Zero-temperature energy and uncertainty for n = 0.8, T = 0, for a range of interaction strengths U, obtained from 
AFQMC (constrained path), DMET, DMRG, FN, MRPHF and UCCSD. Where extrapolations to the TL are not available, 
finite size geometries are listed in lieu of uncertainties. 


n=0.875 

U 

2 

4 

6 

8 

AFQMC 

- 

- 

-1.026 

0.001 

- 

- 

-0.766 

0.001 

DMET 

-1.2721 

0.0006 

-1.031 

0.003 

-0.863 

0.013 

-0.749 

0.007 

DMRG 

- 

- 

-1.028 

[6 X (X)] 

- 

- 

-0.759 

0.004 

FN 

-1.270 

0.002 

-1.0225 

0.0015 

-0.854 

0.002 

-0.749 

0.002 

MRPHF 

-1.2855 

[16 X 4] 

-1.0195 

[16 X 4] 

-0.8318 

[16 X 4] 

-0.7094 

[16 X 4] 

UCCSD 

-1.2667 [16 X 12] 

-1.0093 

[16 X 12] 

-0.8298 [16 X 12] 

-0.7123 [16 X 12] 

UCCSDT* 

-1.2681 

[16 X 4] 

-1.0253 

[16 X 4] 

-0.8570 

[16 X 4] 

-0.7434 

[16 X 4] 


TABLE V. Zero-temperature energy and uncertainty for n = 0.875, T — 0, for a 
range of interaction strengths U, obtained from AFQMC (constrained path), DMET, 
DMRG, FN, MRPHF and UCCSD. Where extrapolations to the TL are not available, 
finite size geometries are listed in lieu of uncertainties. 

Due to strong spatial inhomogeneity observed at this filling, the TL extrapolated 
number in the table is not meaningful, as different cluster shapes show different 
orders. Out of the clusters used here, the 8x2 impurity cluster is probably the 
best estimate with E/t=-0.755 {+/- 0.007). An exhaustive study of the DMET 
cluster size dependence at this filling will be carried out in the future. 
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